linse/sicspsi
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- Updated Makefiles

Browse Source 
- Moved TAS code to psi
- Updated programmers documentation
This commit is contained in:
cvs
2003-06-30 11:51:39 +00:00
parent 9696421d66
commit cd13637987
10 changed files with 3439 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
hardsup/el734_utility.c
View File

@ -2048,6 +2048,9 @@
/* Check the responses carefully. The 3rd response should
** be the device identifier (if to be checked). The 4th
** response should be a hex integer.
** MK: But first we check for *ES, emergency stop.
** The desired behaviour is that we continue in this
** case and hope for the best.
*/
rply_ptr1 = rply_ptr2 = rply_ptr3 = NULL;
rply_ptr0 = AsynSrv_GetReply (
@ -2064,6 +2067,12 @@
if (rply_ptr2 == NULL) rply_ptr2 = "?no_response";
if (rply_ptr3 == NULL) rply_ptr3 = "?no_response";
if(strstr(rply_ptr0,"*ES") != NULL){
*handle = my_handle;
EL734_call_depth--;
return True;
}
if (*rply_ptr1 == '?') rply_ptr0 = rply_ptr1;
if (*rply_ptr2 == '?') rply_ptr0 = rply_ptr2;
if (*rply_ptr3 == '?') rply_ptr0 = rply_ptr3;

2
hardsup/makefile_alpha
View File

@ -10,6 +10,6 @@
#==========================================================================
CC = cc
CFLAGS = -std1 -g $(DFORTIFY) -I$(SRC).. -I$(SRC).
CFLAGS = -std1 -I../.. -g $(DFORTIFY) -I$(SRC).. -I$(SRC).
include make_gen

3
make_gen
View File

@ -13,7 +13,8 @@ OBJ=psi.o buffer.o ruli.o dmc.o nxsans.o nextrics.o sps.o pimotor.o \
ecbcounter.o el737driv.o sinqhmdriv.o tdchm.o velodorn.o \
velodornier.o docho.o sanscook.o tecsdriv.o itc4driv.o itc4.o\
bruker.o ltc11.o A1931.o dilludriv.o eurodriv.o slsmagnet.o \
el755driv.o amorscan.o serial.o scontroller.o
el755driv.o amorscan.o serial.o scontroller.o t_update.o \
t_rlp.o t_conv.o
libpsi.a: $(OBJ)
- rm libpsi.a

17
nextrics.c
View File

@ -562,6 +562,12 @@ name of hkl object holding crystallographic information
SCWrite(pCon,"ERROR: failed to write instrument name",eError);
}
/* Collimator */
SNXSPutMotor(pServ->pSics,pCon,hfil,self->pDict,
"cex1","cex1");
SNXSPutMotor(pServ->pSics,pCon,hfil,self->pDict,
"cex2","cex2");
/* Monochromator */
pVar = NULL;
@ -1029,8 +1035,17 @@ name of hkl object holding crystallographic information
iRet = NXDaliaslink(hfil,self->pDict,"inst0","mono0");
if(iRet != NX_OK)
{
SCWrite(pCon,"WARNING: cannot link against monochromator group",eWarning);
SCWrite(pCon,"WARNING: cannot link against monochromator group",
eWarning);
}
iRet = NXDaliaslink(hfil,self->pDict,"inst0","coll0");
if(iRet != NX_OK)
{
SCWrite(pCon,"WARNING: cannot link against collimator group",
eWarning);
}
/*
links in detector group

838
t_conv.c Normal file
View File

@ -0,0 +1,838 @@
/* t_conv.f -- translated by f2c (version 20000817).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Common Block Declarations */
struct {
integer icrm, icra, iclm;
real cm1rx, cm2rx, ca1rx, ca2rx, rmmin, rmmax, ramin, ramax;
integer inx;
logical at_sinq__;
} curve_;
#define curve_1 curve_
/* Table of constant values */
static integer c__0 = 0;
static integer c__1 = 1;
static doublereal c_b10 = 1.;
static doublereal c_b12 = 360.;
/* ------------------------------------------------------------------------ */
/* slightly edited version for inclusion into SICS */
/* Mark Koennecke, November 2000 */
/* Found that ERRESO looks error messages up in a 2D array. Modified IER */
/* values to refer to a 1D array. */
/* Mark Koennecke, January 2002 */
/* ------------------------------------------------------------------------- */
/* Subroutine */ int inicurve_(integer *midx, real *mrx1, real *mrx2, integer
*aidx, real *arx1, real *arx2, real *mmin, real *mmax, real *amin,
real *amax)
{
/* Initializes a common with the curvature parameters. */
/* In: monochrmoator curvatuure motor index and parameters */
/* In: analyzer curvature motor index + parameters */
curve_1.icrm = *midx;
curve_1.icra = *aidx;
curve_1.cm1rx = *mrx1;
curve_1.cm2rx = *mrx2;
curve_1.ca1rx = *arx1;
curve_1.ca2rx = *arx2;
curve_1.rmmin = *mmin;
curve_1.rmmax = *mmax;
curve_1.ramin = *amin;
curve_1.ramax = *amax;
curve_1.inx = 0;
curve_1.at_sinq__ = TRUE_;
curve_1.iclm = 0;
return 0;
} /* inicurve_ */
/* -------------------------------------------------------------------------- */
/* Subroutine */ int t_conv__(real *ei, real *aki, real *ef, real *akf, real *
qhkl, real *en, real *hx, real *hy, real *hz, integer *if1, integer *
if2, logical *ldk, logical *ldh, logical *ldf, logical *lpa, real *dm,
real *da, real *helm, real *f1h, real *f1v, real *f2h, real *f2v,
real *f, integer *ifx, integer *iss, integer *ism, integer *isa, real
*t_a__, real *t_rm__, real *t_alm__, real *t_ra__, real *qm, logical *
ldra, logical *ldr_rm__, logical *ldr_alm__, logical *ldr_ra__, real *
p_ih__, real *c_ih__, doublereal *a4, integer *ier)
{
/* System generated locals */
doublereal d__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static doublereal edef[2], dakf, daki;
static integer imod;
extern /* Subroutine */ int sam_case__(doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *, integer *);
static integer i__;
static doublereal akdef[2];
extern /* Subroutine */ int helm_case__(real *, real *, real *, real *,
real *, real *, real *, doublereal *, real *, real *, integer *);
static doublereal dqhkl[3];
extern /* Subroutine */ int flip_case__(integer *, integer *, real *,
real *, real *, real *, real *, real *, real *, integer *);
static logical lmoan[2];
static doublereal a1, a2, a3, a5, a6;
static integer id;
static doublereal ra;
extern /* Subroutine */ int rl2spv_(doublereal *, doublereal *,
doublereal *, doublereal *, integer *);
static integer iq;
static doublereal rm;
static logical lqhkle;
extern /* Subroutine */ int erreso_(integer *, integer *);
static doublereal dda, ala, def, dei, ddm, alm, dqm, dqs, dqt[3];
extern /* Subroutine */ int ex_case__(doublereal *, integer *, doublereal
*, doublereal *, doublereal *, doublereal *, doublereal *,
integer *, integer *);
/* ================= */
/* dec$ ident 'V01D' */
/* ----------------------------------------------------------------------- */
/* Other routines in this file: */
/* SUBROUTINE EX_CASE (DX,ISX,AKX,AX1,AX2,RX,ALX,IER) */
/* SUBROUTINE SAM_CASE (QT,QM,QS,AKI,AKF,AX3,AX4,ISS,IER) */
/* SUBROUTINE HELM_CASE (HX,HY,HZ,P_IH,AKI,AKF,A4,QM,HELM,IER) */
/* SUBROUTINE FLIP_CASE (IF1,IF2,P_IH,F1V,F1H,F2V,F2H,AKI,AKF,IER) */
/* ----------------------------------------------------------------------- */
/* INPUTS */
/* EI,AKI,EF,AKF,QHKL,EN,HX,HY,HZ : POTENTIAL TARGETS */
/* IF1,IF2 Status of flippers On (1) Off (0) */
/* LDK(8) LOGICAL INDICATING IF (ENERGY,K OR Q) ARE DRIVEN */
/* LDH,LDF LOGICAL INDICATING IF (HX,HY,HZ) OR (F1,F2) ARE DRIVEN */
/* configuration of the machine */
/* LPA LOGICAL TRUE IF MACHINE IN POLARIZATION MODE */
/* DM,DA,HELM,F1H,F1V,F2H,F2V,F,IFX,ISS,ISM,ISA,QM (F ENERGY UNIT) */
/* OUTPUTs */
/* T_A TARGETS OF ANGLES A1-A6 */
/* T_RM,T_ALM TARGETS OF RM ,LM */
/* T_RA TARGET OF RA */
/* QM TARGETS OF QM */
/* LDRA LOGICAL INDICATING WHICH ANGLES ARE TO BE DRIVEN */
/* LDR_RM LOGICAL INDICATING IF RM (mono curve) IS TO BE DRIVEN */
/* LDR_ALM LOGICAL INDICATING IF ALM (mono transl) IS TO BE DRIVEN */
/* LDR_RA LOGICAL INDICATING IF RA (anal curve) IS TO BE DRIVEN */
/* P_IH TARGETS OF CURRENTS FOR FLIPPERS AND HELMOTZ (8 CURRENTS) */
/* C_IH CONVERSION FACTORS FOR HELMOTZ (4 CURRENTS) */
/* SPECIAL OUTPUTS */
/* TARGET OF EI(EF) IS UPDATED IS KI(KF) IS DRIVEN */
/* TARGET OF VARIABLE ENERGY IS UPDATED IF EN IS DRIVEN */
/* ----------------------------------------------------------------------- */
/* File [MAD.SRC]T_CONV.FOR */
/* include 'curve.inc' */
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
/* ----------------------------------------------------------------------- */
/* LOCAL VARIABLES */
/* ----------------------------------------------------------------------- */
/* SET UP */
/* IMOD INDEX FOR ERROR TREATMENAT BY ERRESO */
/* LDQHKLE : LOGICAL INDICATING THAT WE ARE DEALING WITH A MOVE */
/* IN RECIPROCICAL SPACE */
/* WE REMAP THE ENERGY PB AS FIXED ENERGY IN EDEF(1) */
/* AND VARIABLE ENERGY IN EDEF(2) */
/* IF ISA IS NUL SET IFX TO 1 AND PUT EF,KF, EQUAL TO EI,KI */
/* Parameter adjustments */
--c_ih__;
--p_ih__;
--ldra;
--t_a__;
--ldk;
--qhkl;
/* Function Body */
imod = 3;
ddm = *dm;
dda = *da;
for (i__ = 1; i__ <= 2; ++i__) {
lmoan[i__ - 1] = FALSE_;
}
lqhkle = FALSE_;
for (iq = 5; iq <= 8; ++iq) {
lqhkle = lqhkle || ldk[iq];
}
daki = *aki;
dakf = *akf;
if (*isa == 0) {
*ifx = 1;
}
edef[*ifx - 1] = *ei;
akdef[*ifx - 1] = *aki;
edef[3 - *ifx - 1] = *ef;
akdef[3 - *ifx - 1] = *akf;
if (*isa == 0) {
edef[1] = edef[0];
akdef[1] = akdef[0];
ldk[3] = TRUE_;
ldk[4] = TRUE_;
t_a__[5] = 0.f;
t_a__[6] = 0.f;
ldra[5] = TRUE_;
ldra[6] = TRUE_;
}
/* ----------------------------------------------------------------------- */
/* FIRST TAKE IN ACCOUNT THE FIXED ENERGY PB */
if (ldk[(*ifx << 1) - 1] || ldk[*ifx * 2]) {
lmoan[*ifx - 1] = TRUE_;
if (ldk[(*ifx << 1) - 1]) {
*ier = 9;
if (edef[0] < .1f) {
goto L999;
}
*ier = 0;
akdef[0] = sqrt(edef[0] / *f);
} else {
*ier = 9;
if (akdef[0] < .1f) {
goto L999;
}
*ier = 0;
/* Computing 2nd power */
d__1 = akdef[0];
edef[0] = *f * (d__1 * d__1);
}
}
/* ----------------------------------------------------------------------- */
/* NOW TAKE IN ACCOUNT THE VARIABLE ENERGY PB */
/* VARIABLE ENERGUY IS DRIVEN EITHER EXPLICITLY */
/* E.G. BY DRIVING EI OR KI WITH IFX=2 */
/* ( AND WE MUST CALCULATE EN FROM EVAR) */
/* THE RULE IS : EI=EF+EN : EN IS THE ENERGY LOSS OF NEUTRONS */
/* OR ENERGY GAIN OF SAMPLE */
/* OR IMPLICITLY BY DRIVING THE TRANSFERED ENERGY EN */
/* ( AND WE MUST CALCULATE EVAR FROM EN) */
/* IF KI IS CONSTANT USE THE CURRENT VALUE CONTAINED IN POSN ARRAY */
/* TO CALCULATE THE NON-"CONSTANT" K. */
/* IF KF IS CONSTANT USE ALWAYS THE VALUE IN TARGET AND */
/* DO A DRIVE OF KF TO KEEP A5/A6 IN RIGHT POSITION */
if (ldk[5 - (*ifx << 1)] || ldk[6 - (*ifx << 1)]) {
lmoan[3 - *ifx - 1] = TRUE_;
if (ldk[5 - (*ifx << 1)]) {
*ier = 9;
if (edef[1] < 1e-4f) {
goto L999;
}
*ier = 0;
akdef[1] = sqrt(edef[1] / *f);
} else {
*ier = 9;
if (akdef[1] < 1e-4f) {
goto L999;
}
*ier = 0;
/* Computing 2nd power */
d__1 = akdef[1];
edef[1] = *f * (d__1 * d__1);
}
*en = (3 - (*ifx << 1)) * (edef[0] - edef[1]);
} else if (lqhkle) {
lmoan[3 - *ifx - 1] = TRUE_;
edef[1] = edef[0] + ((*ifx << 1) - 3) * *en;
*ier = 9;
if (edef[1] < 1e-4f) {
goto L999;
}
*ier = 0;
akdef[1] = sqrt(edef[1] / *f);
}
/* ----------------------------------------------------------------------- */
/* CALCULATE MONOCHROMATOR AND ANALYSER ANGLES */
if (lmoan[0]) {
dei = edef[*ifx - 1];
daki = akdef[*ifx - 1];
ex_case__(&ddm, ism, &daki, &a1, &a2, &rm, &alm, &c__0, ier);
if (*ier == 0) {
*aki = daki;
*ei = dei;
t_a__[1] = a1;
t_a__[2] = a2;
ldra[1] = TRUE_;
ldra[2] = TRUE_;
if (curve_1.icrm != 0) {
*t_rm__ = rm;
*ldr_rm__ = TRUE_;
}
if (curve_1.iclm != 0 && curve_1.inx != 0) {
*t_alm__ = alm;
*ldr_alm__ = TRUE_;
}
} else {
*ier += 8;
goto L999;
}
}
if (lmoan[1]) {
def = edef[3 - *ifx - 1];
dakf = akdef[3 - *ifx - 1];
ex_case__(&dda, isa, &dakf, &a5, &a6, &ra, &ala, &c__1, ier);
if (*ier == 0) {
*akf = dakf;
*ef = def;
t_a__[5] = a5;
t_a__[6] = a6;
ldra[5] = TRUE_;
ldra[6] = TRUE_;
if (curve_1.icra != 0) {
*t_ra__ = ra;
*ldr_ra__ = TRUE_;
}
} else {
*ier += 8;
goto L999;
}
}
/* ----------------------------------------------------------------------- */
/* USE (QH,QK,QL) TO CALCULATE A3 AND A4 */
/* CALCULATE Q1 AND Q2 IN SCATTERING PLANE */
imod = 2;
if (lqhkle) {
for (id = 1; id <= 3; ++id) {
dqhkl[id - 1] = qhkl[id];
}
rl2spv_(dqhkl, dqt, &dqm, &dqs, ier);
if (*ier != 0) {
goto L999;
}
sam_case__(dqt, &dqm, &dqs, &daki, &dakf, &a3, a4, iss, ier);
if (*ier == 0) {
t_a__[3] = a3;
t_a__[4] = *a4;
ldra[3] = TRUE_;
ldra[4] = TRUE_;
*qm = dqm;
} else {
*ier += 4;
goto L999;
}
}
/* ----------------------------------------------------------------------- */
/* DEAL WITH FLIPPERS AND HELMOTZ COILS IF LPA */
if (*lpa) {
if (*ldf) {
flip_case__(if1, if2, &p_ih__[1], f1v, f1h, f2v, f2h, aki, akf,
ier);
}
if (*ldh) {
helm_case__(hx, hy, hz, &p_ih__[1], &c_ih__[1], aki, akf, a4, qm,
helm, ier);
}
}
/* ----------------------------------------------------------------------- */
L999:
if (*ier != 0) {
erreso_(&imod, ier);
}
return 0;
} /* t_conv__ */
/* Subroutine */ int ex_case__(doublereal *dx, integer *isx, doublereal *akx,
doublereal *ax1, doublereal *ax2, doublereal *rx, doublereal *alx,
integer *mon_or_anal__, integer *ier)
{
/* System generated locals */
doublereal d__1;
/* Builtin functions */
double asin(doublereal), sin(doublereal), cos(doublereal), sqrt(
doublereal);
/* Local variables */
static integer indx;
static doublereal dcl1r, dc1rx, dc2rx, drmin, drmax, my_rx__, arg;
/* ================== */
/* CALCULATE ANGLES ON MONO/ANALYSER */
/* CALCULATE AX1 AX2 */
/* CALCULATE RX = MONO OR ANAL CURVATURE AND LM = MONO POSIT FOR IN8 */
/* INPUTS */
/* DX D-SPACINGS */
/* ISX SENS OF SCATTERING ON CRYSTAL. If =0, this is probably */
/* a 3-axis instr. in simulated 2-axis mode and the */
/* calculation is for the scattering at the analyser. */
/* In this case, we set AX1 = AX2 = 0 which gives a */
/* "straight-through" setting of A5 & A6 (because of */
/* a simultaneous 90 degree zero offset for A5 -- this */
/* is a bit of a hack, if you ask me!). */
/* AKX TARGET OF MOMENTUM */
/* MON_OR_ANAL =0 if calculation is for mono. */
/* =1 if calculation is for anal. */
/* OUTPUTS */
/* AX1 AX2 THETA 2*THETA ANGLES */
/* RX MONO OR ANALYSER CURVATURE */
/* ALX SPECIAL TRANSLATION FOR IN8 */
/* IER */
/* 1 'KI OR KF CANNOT BE OBTAINED CHECK D-SPACINGS', */
/* 2 'KI OR KF TOO SMALL', */
/* 3 'KI OR KF TOO BIG', */
/* ----------------------------------------------------------------------- */
/* Part of T_CONV.FOR */
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
/* ----------------------------------------------------------------------- */
/* include 'curve.inc' */
/* include 'motdef.inc' */
/* include 'iolsddef.inc' */
/* real*4 tbut(5,NBMOT) */
/* equivalence (rbut, tbut) */
/* ----------------------------------------------------------------------- */
/* LOCAL VAR */
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
*ier = 0;
*ax1 = 0.f;
*ax2 = 0.f;
*rx = 0.f;
*alx = 0.f;
/* ---------------------------------------------------------------------- */
/* Check validity of inputs. */
if (*dx < .1) {
*ier = 1;
}
if (*akx < .1) {
*ier = 2;
}
arg = 3.1415926535897932384626433832795 / (*dx * *akx);
if (abs(arg) > 1.f) {
*ier = 3;
}
if (*ier != 0) {
goto L999;
}
/* ---------------------------------------------------------------------- */
if (*mon_or_anal__ == 0) {
/* Use monochr or anal params? */
indx = curve_1.icrm;
/* Monochr, so set up params. */
dc1rx = curve_1.cm1rx;
dc2rx = curve_1.cm2rx;
dcl1r = (doublereal) curve_1.iclm;
drmin = curve_1.rmmin;
drmax = curve_1.rmmax;
} else {
indx = curve_1.icra;
/* Analyser, so set up params. */
dc1rx = curve_1.ca1rx;
/* There is no ALX in this case. */
dc2rx = curve_1.ca2rx;
drmin = curve_1.ramin;
drmax = curve_1.ramax;
}
/* if (indx .ne. 0) then ! Include zero-offset in min/max */
/* drmin = drmin + tbut(3,indx) */
/* drmax = drmax + tbut(3,indx) */
/* if (drmin .lt. tbut(1,indx)) drmin = tbut(1,indx) */
/* if (drmax .gt. tbut(2,indx)) drmax = tbut(2,indx) */
/* endif */
/* ----------------------------------------------------------------------- */
/* Calculation of the two angles */
if (*isx == 0) {
/* "Straight-through" mode? */
*ax1 = 0.f;
/* Yes. */
*ax2 = 0.f;
*rx = drmin;
*alx = 0.f;
return 0;
}
*ax1 = asin(arg) * *isx * 57.29577951308232087679815481410517;
*ax2 = *ax1 * 2.;
/* ----------------------------------------------------------------------- */
/* Calculation of mono curvature RM or analyser curvature RA */
/* Standard law is: */
/* For monochr: */
/* CM1RX + CM2RX/SIN(abs(A1)/RD) */
/* For analyser: */
/* CA1RX + CA2RX*SIN(abs(A5)/RD) */
/* CM1RX/CM2RX/CA1RX/CA2RX are parameters depending on monochr/analyser and */
/* instrument. They are read from CURVE.INI in routine SETUP_MOT_CURVE. */
/* e.g. cm1rx = .47 */
/* cm2rx = .244 */
/* rmmin = 0. */
/* rmmax = 20. */
/* ----------------------------------------------------------------------- */
if (*mon_or_anal__ == 0) {
/* Monochr or analyser? */
if (curve_1.inx != 0) {
/* Monochr. Is there a translation? */
if (curve_1.iclm != 0) {
/* Yes, IN8 case. If there's a .. */
*alx = dcl1r / sin(*ax2 / 57.29577951308232087679815481410517)
* cos(*ax2 / 57.29577951308232087679815481410517);
/* .. motor, do the .. */
*rx = dc2rx * sqrt(sin(abs(*ax2) /
57.29577951308232087679815481410517)) - dc1rx;
/* .. calculation. */
/* Computing MIN */
d__1 = max(*rx,drmin);
*rx = min(d__1,drmax);
return 0;
}
} else {
/* Not IN8 case so, .. */
my_rx__ = dc1rx + dc2rx / sin(abs(*ax1) /
57.29577951308232087679815481410517);
/* .. simply calculate. */
}
} else {
/* Analyser. */
my_rx__ = dc1rx + dc2rx * sin(abs(*ax1) /
57.29577951308232087679815481410517);
/* Simply calculate. */
}
if (indx != 0) {
/* If there's a motor, return the curvature. */
/* Computing MIN */
d__1 = max(my_rx__,drmin);
*rx = min(d__1,drmax);
/* if (rx .ne. my_rx) then */
/* write (iolun, 101, iostat=ios) motnam(indx), my_rx */
/* 101 format (' Warning -- ', a8, 'curvature restricted by low ', */
/* + 'or high limits!'/ */
/* + ' Calculated curvature was', f9.2) */
/* endif */
}
/* ----------------------------------------------------------------------- */
L999:
return 0;
} /* ex_case__ */
/* Subroutine */ int sam_case__(doublereal *qt, doublereal *qm, doublereal *
qs, doublereal *aki, doublereal *akf, doublereal *ax3, doublereal *
ax4, integer *iss, integer *ier)
{
/* System generated locals */
doublereal d__1, d__2;
/* Builtin functions */
double acos(doublereal), atan2(doublereal, doublereal), d_sign(doublereal
*, doublereal *), d_mod(doublereal *, doublereal *);
/* Local variables */
static doublereal arg, sax3;
/* =================== */
/* DEAL WITH SAMPLE ANGLES CALCULATION FROM Q VECTOR IN C-N PLANE */
/* CALCULATE A3 AND A4 */
/* INPUTS */
/* QT Q-VECTOR IN SCATTERING PLANE */
/* QM,QS Q MODULUS AND QMODULUS SQUARED */
/* AKI,AKF MOMEMTA ON MONO AND ANYLSER */
/* ISS SENS OF SCATTERING ON SAMPLE */
/* OUTPUTS */
/* AX3 AX4 ANGLES ON SAMPLES */
/* IER SAME ERROR AS RL2SPV */
/* IER */
/* 1 'MATRIX S NOT OK', */
/* 2 'Q NOT IN SCATTERING PLANE', */
/* 3 'Q MODULUS TOO SMALL', */
/* 4 'Q MODULUS TOO BIG', */
/* ----------------------------------------------------------------------- */
/* Part of T_CONV.FOR */
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
/* ----------------------------------------------------------------------- */
/* Local variables */
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
/* Parameter adjustments */
--qt;
/* Function Body */
*ier = 0;
if (abs(*qs) < 1e-6 || abs(*qm) < .001) {
*ier = 3;
goto L999;
}
/* ----------------------------------------------------------------------- */
/* CALCULATE A3 AND MOVE IT INTHE -180 ,+180 INTERVAL */
/* Computing 2nd power */
d__1 = *aki;
/* Computing 2nd power */
d__2 = *akf;
arg = (d__1 * d__1 + d__2 * d__2 - *qs) / (*aki * 2. * *akf);
if (abs(arg) > 1.) {
*ier = 4;
goto L999;
} else {
*ax4 = acos(arg) * *iss * 57.29577951308232087679815481410517;
}
/* Computing 2nd power */
d__1 = *akf;
/* Computing 2nd power */
d__2 = *aki;
*ax3 = (-atan2(qt[2], qt[1]) - acos((d__1 * d__1 - *qs - d__2 * d__2) / (*
qm * -2. * *aki)) * d_sign(&c_b10, ax4)) *
57.29577951308232087679815481410517;
sax3 = d_sign(&c_b10, ax3);
d__1 = *ax3 + sax3 * 180.;
*ax3 = d_mod(&d__1, &c_b12) - sax3 * 180.;
/* IF(LPLATE) AX3 = -ATAN(SIN(AX4/RD)/(LSA*TAN(AX5/RD)/(ALMS*C */
/* 1 TAN(AX1/RD))*(AKI/KF)**2-COS(AX4/RD)))*RD !PLATE FOCALIZATION OPTION */
/* IF(AXX3 .GT. 180.D0) AX3 = AX3-360.D0 */
/* IF( A3 .LT. -180.D0) AX3 = AX3+360.D0 */
/* IF(LPLATE .AND. (A3 .GT. 0.0)) AX3 = AX3-180 */
/* C----------------------------------------------------------------------- */
L999:
return 0;
} /* sam_case__ */
/* Subroutine */ int helm_case__(real *hx, real *hy, real *hz, real *t_ih__,
real *c_ih__, real *aki, real *akf, doublereal *a4, real *qm, real *
helm, integer *ier)
{
/* System generated locals */
integer i__1;
real r__1, r__2;
/* Builtin functions */
double cos(doublereal), sin(doublereal), atan2(doublereal, doublereal),
sqrt(doublereal);
/* Local variables */
static doublereal hrad, hdir, qpar, hdir2;
static integer ncoef;
static doublereal qperp;
static integer ic;
static doublereal phi;
static logical at_sinq__;
/* ==================== */
/* DEAL WITH HELMOTZ COIL FIELD CALCULATIONS */
/* CALCULATE HX HY HZ */
/* ----------------------------------------------------------------------- */
/* 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. */
/* ----------------------------------------------------------------------- */
/* Part of T_CONV.FOR */
/* include 'common_sinq.inc' */
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
/* ----------------------------------------------------------------------- */
/* Local variables */
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
/* Parameter adjustments */
--c_ih__;
--t_ih__;
/* Function Body */
at_sinq__ = TRUE_;
ncoef = 4;
if (at_sinq__) {
ncoef = 3;
}
*ier = 1;
if (dabs(*qm) < 1e-4) {
goto L999;
}
*ier = 0;
i__1 = ncoef;
for (ic = 1; ic <= i__1; ++ic) {
if (c_ih__[ic] < 1e-4) {
*ier = 2;
}
}
if (*ier != 0) {
goto L999;
}
/* ----------------------------------------------------------------------- */
/* CALCULATE MODULE AND ANGLES OF IN PLANE FIELD H */
/* PHI ! Angle between Q and KI (in radians) */
/* HRAD ! Radial comp. of H */
/* HDIR ! Direction of H relative to PHI (in radians) */
/* HDIR2 ! Angle between field and axis of Coil 1 (in radians) */
qpar = *aki - *akf * cos(*a4 / 57.29577951308232087679815481410517);
qperp = -(*akf) * sin(*a4 / 57.29577951308232087679815481410517);
if (abs(qpar) > .001 && abs(qperp) > .001) {
phi = atan2((abs(qperp)), (abs(qpar)));
if (qpar > 0. && qperp < 0.) {
phi = -phi;
} else if (qpar < 0. && qperp > 0.) {
phi = 3.1415926535897932384626433832795 - phi;
} else if (qpar < 0. && qperp < 0.) {
phi += -3.1415926535897932384626433832795;
}
} else if (abs(qpar) > .001) {
if (qpar >= 0.f) {
phi = 0.f;
}
if (qpar < 0.f) {
phi = 3.1415926535897932384626433832795;
}
} else if (abs(qperp) > .001) {
if (qperp >= 0.f) {
phi = 1.5707963267948966;
}
if (qperp < 0.f) {
phi = -1.5707963267948966;
}
} else {
phi = 0.f;
}
/* Computing 2nd power */
r__1 = *hx;
/* Computing 2nd power */
r__2 = *hy;
hrad = sqrt(r__1 * r__1 + r__2 * r__2);
if (dabs(*hx) > .001 && dabs(*hy) > .001) {
hdir = atan2((dabs(*hy)), (dabs(*hx)));
if (*hx > 0.f && *hy < 0.f) {
hdir = -hdir;
} else if (*hx < 0.f && *hy > 0.f) {
hdir = 3.1415926535897932384626433832795 - hdir;
} else if (*hx < 0.f && *hy < 0.f) {
hdir += -3.1415926535897932384626433832795;
}
} else if (dabs(*hx) > .001) {
if (*hx >= 0.f) {
hdir = 0.f;
}
if (*hx < 0.f) {
hdir = 3.1415926535897932384626433832795;
}
} else if (dabs(*hy) > .001) {
if (*hy >= 0.f) {
hdir = 1.5707963267948966;
}
if (*hy < 0.f) {
hdir = -1.5707963267948966;
}
} else {
hdir = 0.f;
}
hdir2 = hdir + phi - *helm / 57.29577951308232087679815481410517;
/* ----------------------------------------------------------------------- */
/* !CALC CURRENTS */
/* !POSITION OF PSP FOR COIL I */
if (! at_sinq__) {
hdir2 += 1.5707963267948966;
/* ??? */
for (ic = 1; ic <= 3; ++ic) {
t_ih__[ic + 4] = cos(hdir2 + (ic - 1) * 2.f *
3.1415926535897932384626433832795 / 3.f) * hrad / c_ih__[
ic] / 1.5f;
}
t_ih__[8] = *hz / c_ih__[4];
} else {
t_ih__[5] = cos(hdir2) * hrad / c_ih__[1];
t_ih__[6] = sin(hdir2) * hrad / c_ih__[2];
t_ih__[7] = *hz / c_ih__[3];
}
/* ----------------------------------------------------------------------- */
L999:
return 0;
} /* helm_case__ */
/* Subroutine */ int flip_case__(integer *if1, integer *if2, real *t_ih__,
real *f1v, real *f1h, real *f2v, real *f2h, real *aki, real *akf,
integer *ier)
{
/* ==================== */
/* DEAL WITH FLIPPER COIL CALCULATIONS */
/* CALCULATE P_IF CURRENTS FOR THE TWO FLIPPERS */
/* ----------------------------------------------------------------------- */
/* Define the dummy arguments */
/* Part of T_CONV.FOR */
/* ----------------------------------------------------------------------- */
/* INIT AND TEST */
/* Parameter adjustments */
--t_ih__;
/* Function Body */
*ier = 0;
/* ----------------------------------------------------------------------- */
if (*if1 == 1) {
t_ih__[1] = *f1v;
t_ih__[2] = *aki * *f1h;
} else {
t_ih__[1] = 0.f;
t_ih__[2] = 0.f;
}
if (*if2 == 1) {
t_ih__[3] = *f2v;
t_ih__[4] = *akf * *f2h;
} else {
t_ih__[3] = 0.f;
t_ih__[4] = 0.f;
}
/* ----------------------------------------------------------------------- */
/* L999: */
return 0;
} /* flip_case__ */

690
t_conv.f Executable file
View File

@ -0,0 +1,690 @@
C------------------------------------------------------------------------
C slightly edited version for inclusion into SICS
C
C Mark Koennecke, November 2000
C
C Found that ERRESO looks error messages up in a 2D array. Modified IER
C values to refer to a 1D array.
C
C Mark Koennecke, January 2002
C-------------------------------------------------------------------------
SUBROUTINE INICURVE(MIDX, MRX1, MRX2, AIDX, ARX1, ARX2,
+ MMIN, MMAX, AMIN, AMAX)
C
C Initializes a common with the curvature parameters.
C In: monochrmoator curvatuure motor index and parameters
C In: analyzer curvature motor index + parameters
INTEGER MIDX, AIDX
REAL*4 MRX1, MRX2, ARX1, ARX2, MMIN, MMAX, AMIN, AMAX
REAL*4 CM1RX, CM2RX, CA1RX, CA2RX, RMMIN, RMMAX
REAL*4 RAMIN, RAMAX
INTEGER ICRM, ICRA, ICLM, INX
LOGICAL AT_SINQ
COMMON/CURVE/ICRM,ICRA, ICLM, CM1RX, CM2RX, CA1RX, CA2RX,
+ RMMIN, RMMAX, RAMIN, RAMAX, INX, AT_SINQ
ICRM = MIDX
ICRA = AIDX
CM1RX = MRX1
CM2RX = MRX2
CA1RX = ARX1
CA2RX = ARX2
RMMIN = MMIN
RMMAX = MMAX
RAMIN = AMIN
RAMAX = AMAX
INX = 0
AT_SINQ = .TRUE.
ICLM = 0
RETURN
END
C--------------------------------------------------------------------------
SUBROUTINE T_CONV ( ! File [MAD.SRC]T_CONV.FOR
c =================
+ EI, AKI, EF, AKF, QHKL, EN,
+ HX, HY, HZ,
+ IF1, IF2,
+ LDK, LDH, LDF, LPA,
+ DM, DA,
+ HELM, F1H, F1V, F2H, F2V, F,
+ IFX, ISS, ISM, ISA,
+ T_A, T_RM, T_ALM, T_RA, QM,
+ LDRA, LDR_RM, LDR_ALM, LDR_RA,
+ P_IH, C_IH, A4,
+ IER)
c
cdec$ ident 'V01D'
C-----------------------------------------------------------------------
C Other routines in this file:
c
C SUBROUTINE EX_CASE (DX,ISX,AKX,AX1,AX2,RX,ALX,IER)
C SUBROUTINE SAM_CASE (QT,QM,QS,AKI,AKF,AX3,AX4,ISS,IER)
C SUBROUTINE HELM_CASE (HX,HY,HZ,P_IH,AKI,AKF,A4,QM,HELM,IER)
C SUBROUTINE FLIP_CASE (IF1,IF2,P_IH,F1V,F1H,F2V,F2H,AKI,AKF,IER)
C-----------------------------------------------------------------------
C INPUTS
C EI,AKI,EF,AKF,QHKL,EN,HX,HY,HZ : POTENTIAL TARGETS
C IF1,IF2 Status of flippers On (1) Off (0)
C LDK(8) LOGICAL INDICATING IF (ENERGY,K OR Q) ARE DRIVEN
C LDH,LDF LOGICAL INDICATING IF (HX,HY,HZ) OR (F1,F2) ARE DRIVEN
C
C configuration of the machine
c
C LPA LOGICAL TRUE IF MACHINE IN POLARIZATION MODE
C DM,DA,HELM,F1H,F1V,F2H,F2V,F,IFX,ISS,ISM,ISA,QM (F ENERGY UNIT)
C
C OUTPUTs
C T_A TARGETS OF ANGLES A1-A6
C T_RM,T_ALM TARGETS OF RM ,LM
C T_RA TARGET OF RA
C QM TARGETS OF QM
C LDRA LOGICAL INDICATING WHICH ANGLES ARE TO BE DRIVEN
C LDR_RM LOGICAL INDICATING IF RM (mono curve) IS TO BE DRIVEN
C LDR_ALM LOGICAL INDICATING IF ALM (mono transl) IS TO BE DRIVEN
C LDR_RA LOGICAL INDICATING IF RA (anal curve) IS TO BE DRIVEN
C P_IH TARGETS OF CURRENTS FOR FLIPPERS AND HELMOTZ (8 CURRENTS)
C C_IH CONVERSION FACTORS FOR HELMOTZ (4 CURRENTS)
C
C SPECIAL OUTPUTS
C TARGET OF EI(EF) IS UPDATED IS KI(KF) IS DRIVEN
C TARGET OF VARIABLE ENERGY IS UPDATED IF EN IS DRIVEN
C-----------------------------------------------------------------------
implicit none
c
REAL*4 EPS1, EPS4
parameter (EPS1 = 1.D-1)
parameter (EPS4 = 1.D-4)
INTEGER ICRM, ICRA, ICLM, INX
LOGICAL AT_SINQ
REAL*4 CM1RX, CM2RX, CA1RX, CA2RX, RMMIN, RMMAX
REAL*4 RAMIN, RAMAX
COMMON/CURVE/ICRM,ICRA, ICLM, CM1RX, CM2RX, CA1RX, CA2RX,
+ RMMIN, RMMAX, RAMIN, RAMAX, INX, AT_SINQ
c
C include 'curve.inc'
C-----------------------------------------------------------------------
C Define the dummy arguments
c
real*4 ei, aki, ef, akf, qhkl(3), en
real*4 hx, hy, hz
integer*4 if1, if2
logical*4 ldk(8), ldh, ldf, lpa
real*4 dm, da
real*4 helm, f1h, f1v, f2h, f2v, f
integer*4 ifx, iss, ism, isa
real*4 t_a(6), t_rm, t_alm, t_ra, qm
logical*4 ldra(6), ldr_rm, ldr_alm, ldr_ra
real*4 p_ih(8), c_ih(4)
integer*4 ier
C-----------------------------------------------------------------------
C LOCAL VARIABLES
C
integer*4 i, id, imod, iq
logical*4 lmoan(2), lqhkle
c
double precision a1, a2, a3, a4, a5, a6
double precision ala, alm, dakf, daki, dqm, dqs
double precision def, dei
double precision ra, rm
double precision edef(2), akdef(2), dqhkl(3), dqt(3)
double precision ddm, dda
C-----------------------------------------------------------------------
C SET UP
C IMOD INDEX FOR ERROR TREATMENAT BY ERRESO
C LDQHKLE : LOGICAL INDICATING THAT WE ARE DEALING WITH A MOVE
C IN RECIPROCICAL SPACE
C WE REMAP THE ENERGY PB AS FIXED ENERGY IN EDEF(1)
C AND VARIABLE ENERGY IN EDEF(2)
C IF ISA IS NUL SET IFX TO 1 AND PUT EF,KF, EQUAL TO EI,KI
C
IMOD = 3
DDM = DM
DDA = DA
DO I = 1,2
LMOAN(I) = .FALSE.
ENDDO
LQHKLE = .FALSE.
DO IQ = 5,8
LQHKLE = (LQHKLE .OR. LDK(IQ))
ENDDO
DAKI = AKI
DAKF = AKF
IF (ISA .EQ. 0) IFX = 1
EDEF(IFX) = EI
AKDEF(IFX) = AKI
EDEF(3-IFX) = EF
AKDEF(3-IFX) = AKF
IF( ISA .EQ. 0) THEN
EDEF(2) = EDEF(1)
AKDEF(2) = AKDEF(1)
LDK(3) = .TRUE.
LDK(4) = .TRUE.
T_A(5) = 0.
T_A(6) = 0.
LDRA(5) = .TRUE.
LDRA(6) = .TRUE.
ENDIF
C-----------------------------------------------------------------------
C FIRST TAKE IN ACCOUNT THE FIXED ENERGY PB
C
IF (LDK(2*IFX-1) .OR. LDK(2*IFX)) THEN
LMOAN(IFX) = .TRUE.
IF (LDK(2*IFX-1)) THEN
IER = 1 + 8
IF(EDEF(1) .LT. EPS1) GOTO 999
IER = 0
AKDEF(1) = SQRT(EDEF(1)/F)
ELSE
IER = 1 + 8
IF(AKDEF(1) .LT. EPS1) GOTO 999
IER = 0
EDEF(1) = F*AKDEF(1)**2
ENDIF
ENDIF
C-----------------------------------------------------------------------
C NOW TAKE IN ACCOUNT THE VARIABLE ENERGY PB
C VARIABLE ENERGUY IS DRIVEN EITHER EXPLICITLY
C E.G. BY DRIVING EI OR KI WITH IFX=2
C ( AND WE MUST CALCULATE EN FROM EVAR)
C THE RULE IS : EI=EF+EN : EN IS THE ENERGY LOSS OF NEUTRONS
C OR ENERGY GAIN OF SAMPLE
C OR IMPLICITLY BY DRIVING THE TRANSFERED ENERGY EN
C ( AND WE MUST CALCULATE EVAR FROM EN)
C IF KI IS CONSTANT USE THE CURRENT VALUE CONTAINED IN POSN ARRAY
C TO CALCULATE THE NON-"CONSTANT" K.
C IF KF IS CONSTANT USE ALWAYS THE VALUE IN TARGET AND
C DO A DRIVE OF KF TO KEEP A5/A6 IN RIGHT POSITION
C
IF (LDK(5-2*IFX) .OR. LDK(6-2*IFX)) THEN
LMOAN(3-IFX) = .TRUE.
IF (LDK(5-2*IFX)) THEN
IER = 1 + 8
IF(EDEF(2) .LT. EPS4) GOTO 999
IER = 0
AKDEF(2) = SQRT(EDEF(2)/F)
ELSE
IER = 1 + 8
IF(AKDEF(2) .LT. EPS4) GOTO 999
IER = 0
EDEF(2) = F*AKDEF(2)**2
ENDIF
EN = (3-2*IFX)*(EDEF(1)-EDEF(2))
ELSEIF (LQHKLE) THEN
LMOAN(3-IFX) = .TRUE.
EDEF(2) = EDEF(1)+(2*IFX-3)*EN
IER = 1 + 8
IF(EDEF(2) .LT. EPS4) GOTO 999
IER = 0
AKDEF(2) = SQRT(EDEF(2)/F)
ENDIF
C-----------------------------------------------------------------------
C CALCULATE MONOCHROMATOR AND ANALYSER ANGLES
C
IF(LMOAN(1)) THEN
DEI = EDEF(IFX)
DAKI = AKDEF(IFX)
CALL EX_CASE(DDM,ISM,DAKI,A1,A2,RM,ALM,0,IER)
IF (IER .EQ. 0) THEN
AKI = DAKI
EI = DEI
T_A(1) = A1
T_A(2) = A2
LDRA(1) = .TRUE.
LDRA(2) = .TRUE.
if (icrm .ne. 0) then
T_RM = RM
LDR_RM = .TRUE.
endif
if ((iclm .ne. 0) .and. (inx .ne. 0)) then
T_ALM = ALM
LDR_ALM = .TRUE.
endif
ELSE
IER = IER + 8
GOTO 999
ENDIF
ENDIF
IF(LMOAN(2)) THEN
DEF = EDEF(3-IFX)
DAKF = AKDEF(3-IFX)
CALL EX_CASE(DDA,ISA,DAKF,A5,A6,RA,ALA,1,IER)
IF (IER .EQ. 0) THEN
AKF = DAKF
EF = DEF
T_A(5) = A5
T_A(6) = A6
LDRA(5) = .TRUE.
LDRA(6) = .TRUE.
if (icra .ne. 0) then
T_RA = RA
LDR_RA = .TRUE.
endif
ELSE
IER = IER + 8
GOTO 999
ENDIF
ENDIF
C-----------------------------------------------------------------------
C USE (QH,QK,QL) TO CALCULATE A3 AND A4
C CALCULATE Q1 AND Q2 IN SCATTERING PLANE
C
IMOD = 2
IF (LQHKLE) THEN
DO ID = 1,3
DQHKL(ID) = QHKL(ID)
ENDDO
CALL RL2SPV(DQHKL,DQT,DQM,DQS,IER)
IF (IER .NE. 0) GOTO 999
CALL SAM_CASE(DQT,DQM,DQS,DAKI,DAKF,A3,A4,ISS,IER)
IF (IER .EQ. 0) THEN
T_A(3) = A3
T_A(4) = A4
LDRA(3) = .TRUE.
LDRA(4) = .TRUE.
QM = DQM
ELSE
IER = IER + 4
GOTO 999
ENDIF
ENDIF
C-----------------------------------------------------------------------
C DEAL WITH FLIPPERS AND HELMOTZ COILS IF LPA
C
IF (LPA) THEN
IF (LDF) CALL FLIP_CASE(IF1,IF2,P_IH,F1V,F1H,F2V,F2H,
+ AKI,AKF,IER)
IF (LDH) CALL HELM_CASE(HX,HY,HZ,P_IH,C_IH,AKI,AKF,
+ A4,QM,HELM,IER)
endif
C-----------------------------------------------------------------------
999 CONTINUE
IF (IER .NE. 0) CALL ERRESO(IMOD,IER)
RETURN
END
c
SUBROUTINE EX_CASE ( ! Part of T_CONV.FOR
c ==================
+ DX,ISX,AKX,AX1,AX2,RX,ALX,MON_OR_ANAL,IER)
C
C CALCULATE ANGLES ON MONO/ANALYSER
C CALCULATE AX1 AX2
C CALCULATE RX = MONO OR ANAL CURVATURE AND LM = MONO POSIT FOR IN8
C
C INPUTS
C DX D-SPACINGS
C ISX SENS OF SCATTERING ON CRYSTAL. If =0, this is probably
c a 3-axis instr. in simulated 2-axis mode and the
c calculation is for the scattering at the analyser.
c In this case, we set AX1 = AX2 = 0 which gives a
c "straight-through" setting of A5 & A6 (because of
c a simultaneous 90 degree zero offset for A5 -- this
c is a bit of a hack, if you ask me!).
C AKX TARGET OF MOMENTUM
c MON_OR_ANAL =0 if calculation is for mono.
c =1 if calculation is for anal.
C OUTPUTS
C AX1 AX2 THETA 2*THETA ANGLES
C RX MONO OR ANALYSER CURVATURE
C ALX SPECIAL TRANSLATION FOR IN8
C IER
C 1 'KI OR KF CANNOT BE OBTAINED CHECK D-SPACINGS',
C 2 'KI OR KF TOO SMALL',
C 3 'KI OR KF TOO BIG',
C-----------------------------------------------------------------------
implicit none
c
double precision PI, RD, EPS1
PARAMETER (PI = 3.14159265358979323846264338327950D0)
PARAMETER (RD = 57.29577951308232087679815481410517D0)
PARAMETER (EPS1 = 1.D-1)
C-----------------------------------------------------------------------
C Define the dummy arguments
double precision dx
integer*4 isx
double precision akx, ax1, ax2, rx, alx
integer*4 mon_or_anal, ier
C-----------------------------------------------------------------------
C include 'curve.inc'
C include 'motdef.inc'
C include 'iolsddef.inc'
c
INTEGER ICRM, ICRA, ICLM, INX
LOGICAL AT_SINQ
REAL*4 CM1RX, CM2RX, CA1RX, CA2RX, RMMIN, RMMAX
REAL*4 RAMIN, RAMAX
COMMON/CURVE/ICRM,ICRA, ICLM, CM1RX, CM2RX, CA1RX, CA2RX,
+ RMMIN, RMMAX, RAMIN, RAMAX, INX, AT_SINQ
C real*4 tbut(5,NBMOT)
C equivalence (rbut, tbut)
C-----------------------------------------------------------------------
C LOCAL VAR
c
double precision arg, dc1rx, dc2rx, drmin, drmax, dcl1r, my_rx
integer*4 ios, indx
C-----------------------------------------------------------------------
C INIT AND TEST
C
ier = 0
ax1 = 0.0
ax2 = 0.0
rx = 0.0
alx = 0.0
c----------------------------------------------------------------------
c Check validity of inputs.
if (dx .lt. EPS1) ier = 1
if (akx .lt. EPS1) ier = 2
arg = PI/(dx * akx)
if (abs(arg) .gt. 1.0) ier = 3
if (ier .ne. 0) goto 999
c----------------------------------------------------------------------
if (mon_or_anal .eq. 0) then ! Use monochr or anal params?
indx = icrm ! Monochr, so set up params.
dc1rx = cm1rx
dc2rx = cm2rx
dcl1r = ICLM
drmin = rmmin
drmax = rmmax
else
indx = icra ! Analyser, so set up params.
dc1rx = ca1rx ! There is no ALX in this case.
dc2rx = ca2rx
drmin = ramin
drmax = ramax
endif
c
C if (indx .ne. 0) then ! Include zero-offset in min/max
C drmin = drmin + tbut(3,indx)
C drmax = drmax + tbut(3,indx)
C if (drmin .lt. tbut(1,indx)) drmin = tbut(1,indx)
C if (drmax .gt. tbut(2,indx)) drmax = tbut(2,indx)
C endif
c-----------------------------------------------------------------------
c Calculation of the two angles
c
if (isx .eq. 0) then ! "Straight-through" mode?
ax1 = 0.0 ! Yes.
ax2 = 0.0
rx = drmin
alx = 0.0
return
endif
c
ax1 = asin (arg) * isx * rd
ax2 = 2.0d0 * ax1
c-----------------------------------------------------------------------
c Calculation of mono curvature RM or analyser curvature RA
c Standard law is:
c
c For monochr:
c CM1RX + CM2RX/SIN(abs(A1)/RD)
c
c For analyser:
c CA1RX + CA2RX*SIN(abs(A5)/RD)
c
c CM1RX/CM2RX/CA1RX/CA2RX are parameters depending on monochr/analyser and
c instrument. They are read from CURVE.INI in routine SETUP_MOT_CURVE.
c e.g. cm1rx = .47
c cm2rx = .244
c rmmin = 0.
c rmmax = 20.
c-----------------------------------------------------------------------
if (mon_or_anal .eq. 0) then ! Monochr or analyser?
if (inx .ne. 0) then ! Monochr. Is there a translation?
if (iclm .ne. 0) then ! Yes, IN8 case. If there's a ..
alx = (dcl1r/sin(ax2/rd)) * cos(ax2/rd) ! .. motor, do the ..
rx = dc2rx * sqrt(sin(abs(ax2)/rd)) - dc1rx ! .. calculation.
rx = dmin1 (dmax1 (rx, drmin), drmax)
return
endif
else ! Not IN8 case so, ..
my_rx = dc1rx + dc2rx/sin(abs(ax1)/rd) ! .. simply calculate.
endif
else ! Analyser.
my_rx = dc1rx + dc2rx * sin(abs(ax1)/rd) ! Simply calculate.
endif
c
if (indx .ne. 0) then ! If there's a motor, return the curvature.
rx = dmin1 (dmax1 (my_rx, drmin), drmax)
C if (rx .ne. my_rx) then
C write (iolun, 101, iostat=ios) motnam(indx), my_rx
C 101 format (' Warning -- ', a8, 'curvature restricted by low ',
C + 'or high limits!'/
C + ' Calculated curvature was', f9.2)
C endif
endif
c-----------------------------------------------------------------------
999 continue
return
end
c
SUBROUTINE SAM_CASE ( ! Part of T_CONV.FOR
c ===================
+ QT,QM,QS,AKI,AKF,AX3,AX4,ISS,IER)
C
C DEAL WITH SAMPLE ANGLES CALCULATION FROM Q VECTOR IN C-N PLANE
C CALCULATE A3 AND A4
C INPUTS
C QT Q-VECTOR IN SCATTERING PLANE
C QM,QS Q MODULUS AND QMODULUS SQUARED
C AKI,AKF MOMEMTA ON MONO AND ANYLSER
C ISS SENS OF SCATTERING ON SAMPLE
C
C OUTPUTS
C AX3 AX4 ANGLES ON SAMPLES
C IER SAME ERROR AS RL2SPV
C IER
C 1 'MATRIX S NOT OK',
C 2 'Q NOT IN SCATTERING PLANE',
C 3 'Q MODULUS TOO SMALL',
C 4 'Q MODULUS TOO BIG',
C-----------------------------------------------------------------------
implicit none
c
double precision RD, EPS3, EPS6
PARAMETER (RD = 57.29577951308232087679815481410517D0)
PARAMETER (EPS3 = 1.D-3)
PARAMETER (EPS6 = 1.D-6)
C-----------------------------------------------------------------------
C Define the dummy arguments
double precision qt(3)
double precision qm, qs, aki, akf, ax3, ax4
integer*4 iss, ier
C-----------------------------------------------------------------------
c Local variables
double precision arg, sax3
C-----------------------------------------------------------------------
C INIT AND TEST
C
IER = 0
IF ((ABS(QS) .LT. EPS6) .OR. (ABS(QM) .LT. EPS3)) THEN
IER = 3
GOTO 999
ENDIF
C-----------------------------------------------------------------------
C CALCULATE A3 AND MOVE IT INTHE -180 ,+180 INTERVAL
C
ARG = (AKI**2 + AKF**2 - QS)/(2.D0*AKI*AKF)
IF(ABS(ARG) .GT. 1.D0)THEN
IER = 4
GOTO 999
ELSE
AX4 = ACOS(ARG)*ISS*RD
ENDIF
AX3 = (-ATAN2(QT(2),QT(1))-ACOS((AKF**2-QS-AKI**2)/
+ (-2.D0*QM*AKI))*DSIGN(1.D0,AX4))*RD
sax3 = Dsign(1.D0,ax3)
AX3 = DMOD(AX3+sax3*180.D0,360.D0)-sax3*180.D0
C
C IF(LPLATE) AX3 = -ATAN(SIN(AX4/RD)/(LSA*TAN(AX5/RD)/(ALMS*C
C 1 TAN(AX1/RD))*(AKI/KF)**2-COS(AX4/RD)))*RD !PLATE FOCALIZATION OPTION
C IF(AXX3 .GT. 180.D0) AX3 = AX3-360.D0
C IF( A3 .LT. -180.D0) AX3 = AX3+360.D0
C IF(LPLATE .AND. (A3 .GT. 0.0)) AX3 = AX3-180
CC-----------------------------------------------------------------------
999 CONTINUE
RETURN
END
c
SUBROUTINE HELM_CASE ( ! Part of T_CONV.FOR
c ====================
+ HX,HY,HZ,T_IH,C_IH,AKI,AKF,A4,QM,HELM,IER)
C
C DEAL WITH HELMOTZ COIL FIELD CALCULATIONS
C CALCULATE HX HY HZ
C-----------------------------------------------------------------------
c At ILL:
c There are 3 coils for Hx/Hy at 120 degrees to each other.
c
c There is a 4th coil for Hz.
c
c At SINQ:
c There is an Hx coil and an Hy coil (actually each is 4 coils powered
c in series). They are mounted on a ring (SRO). The value of HELM is
c the angle between the Hx coil and ki.
c
c There is a 3rd coil for Hz.
C-----------------------------------------------------------------------
implicit none
c
C include 'common_sinq.inc'
c
double precision PI, RD, EPS3, EPS4
PARAMETER (PI = 3.14159265358979323846264338327950D0)
PARAMETER (RD = 57.29577951308232087679815481410517D0)
PARAMETER (EPS3 = 1.0D-3)
PARAMETER (EPS4 = 1.0D-4)
C-----------------------------------------------------------------------
C Define the dummy arguments
real*4 hx, hy, hz
real*4 t_ih(8)
real*4 c_ih(4)
real*4 aki, akf
double precision a4
real*4 qm, helm
integer*4 ier
LOGICAL AT_SINQ
C-----------------------------------------------------------------------
c Local variables
integer*4 ic, ncoef
double precision hdir, hdir2, hrad, phi, qpar, qperp
C-----------------------------------------------------------------------
C INIT AND TEST
C
AT_SINQ = .TRUE.
ncoef = 4
if (at_sinq) ncoef = 3
c
IER = 1
IF (ABS(QM) .LT. EPS4) goto 999
IER = 0
DO IC = 1,ncoef
IF (C_IH(IC) .LT. EPS4) IER = 2
ENDDO
IF (IER .NE. 0) GOTO 999
C-----------------------------------------------------------------------
C CALCULATE MODULE AND ANGLES OF IN PLANE FIELD H
C PHI ! Angle between Q and KI (in radians)
C HRAD ! Radial comp. of H
C HDIR ! Direction of H relative to PHI (in radians)
C HDIR2 ! Angle between field and axis of Coil 1 (in radians)
C
qpar = aki - akf * cos(a4/RD)
qperp = -akf * sin(a4/RD)
if (abs(qpar) .gt. EPS3 .and. abs(qperp) .gt. EPS3) then
phi = atan2 (abs(qperp), abs(qpar))
if (qpar .gt. 0 .and. qperp .lt. 0) then
phi = -phi
elseif (qpar .lt. 0 .and. qperp .gt. 0) then
phi = PI - phi
elseif (qpar .lt. 0 .and. qperp .lt. 0) then
phi = phi - PI
endif
elseif (abs(qpar) .gt. EPS3) then
if (qpar .ge. 0.0) phi = 0.0
if (qpar .lt. 0.0) phi = PI
elseif (abs(qperp) .gt. EPS3) then
if (qperp .ge. 0.0) phi = 0.5 * PI
if (qperp .lt. 0.0) phi = -0.5 * PI
else
phi = 0.0
endif
c
hrad = sqrt (hx**2 + hy**2)
if (abs(hx) .gt. EPS3 .and. abs(hy) .gt. EPS3) then
hdir = atan2 (abs(hy), abs(hx))
if (hx .gt. 0 .and. hy .lt. 0) then
hdir = -hdir
elseif (hx .lt. 0 .and. hy .gt. 0) then
hdir = PI - hdir
elseif (hx .lt. 0 .and. hy .lt. 0) then
hdir = hdir - PI
endif
elseif (abs(hx) .gt. EPS3) then
if (hx .ge. 0.0) hdir = 0.0
if (hx .lt. 0.0) hdir = PI
elseif (abs(hy) .gt. EPS3) then
if (hy .ge. 0.0) hdir = 0.5 * PI
if (hy .lt. 0.0) hdir = -0.5 * PI
else
hdir = 0.0
endif
c
hdir2 = hdir + phi - (helm/RD)
C-----------------------------------------------------------------------
C !CALC CURRENTS
C !POSITION OF PSP FOR COIL I
C
if (.not. at_sinq) then
hdir2 = hdir2 + 0.5 * PI ! ???
do ic = 1,3
t_ih(ic+4) = cos(hdir2+(ic-1)*2.*PI/3.)*hrad/c_ih(ic)/1.5
enddo
t_ih(8) = hz/c_ih(4)
else
t_ih(5) = cos(hdir2) * hrad/c_ih(1)
t_ih(6) = sin(hdir2) * hrad/c_ih(2)
t_ih(7) = hz/c_ih(3)
endif
C-----------------------------------------------------------------------
999 CONTINUE
RETURN
END
c
SUBROUTINE FLIP_CASE ( ! Part of T_CONV.FOR
C ====================
+ IF1,IF2,T_IH,F1V,F1H,F2V,F2H,AKI,AKF,IER)
C
C DEAL WITH FLIPPER COIL CALCULATIONS
C CALCULATE P_IF CURRENTS FOR THE TWO FLIPPERS
C-----------------------------------------------------------------------
C Define the dummy arguments
integer*4 if1, if2
real*4 t_ih(8)
real*4 f1v, f1h, f2v, f2h
real*4 aki, akf
integer*4 ier
C-----------------------------------------------------------------------
C INIT AND TEST
C
IER = 0
C-----------------------------------------------------------------------
C
IF (IF1 .EQ. 1) THEN
T_IH(1) = F1V
T_IH(2) = AKI*F1H
ELSE
T_IH(1) = 0.
T_IH(2) = 0.
ENDIF
IF (IF2 .EQ. 1) THEN
T_IH(3) = F2V
T_IH(4) = AKF*F2H
ELSE
T_IH(3) = 0.
T_IH(4) = 0.
ENDIF
C-----------------------------------------------------------------------
999 CONTINUE
RETURN
END

516
t_rlp.c Normal file
View File

@ -0,0 +1,516 @@
/* t_rlp.f -- translated by f2c (version 20000817).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Common Block Declarations */
struct {
doublereal s[16] /* was [4][4] */, sinv[16] /* was [4][4] */;
integer iok;
} osolem_;
#define osolem_1 osolem_
/* Subroutine */ int t_rlp__(void)
{
/* ================ */
/* dec$ Ident 'V01A' */
/* ------------------------------------------------------------------ */
/* Updates: */
/* V01A 7-May-1996 DM. Put error output to IOLUN, use IMPLICIT NONE and */
/* get the code indented so that it is readable! */
/* ------------------------------------------------------------------ */
/* Routines to deal with the reciprocical lattice PB */
/* ------------------------------------------------------------------ */
/* Entry points in this file: */
/* SETRLP : CALCULATION OF S AND INVS , ORIENTATION MATRIX */
/* RL2SPV : TRANSFO FROM RECIP LAT TO SCAT PLANE */
/* SP2RLV : TRANSFO FROM SCAT PLANE TO RECIP LAT */
/* INVS : INVERT MATRIX S, GENERATED BY SETRLP. */
/* ERRESO : DEAL ITH ERROR MESSAGES FOR ALL MODULES */
/* SUBROUTINE SETRLP(SAM,IER) */
/* SUBROUTINE RL2SPV(QHKL,QT,QM,QS,IER) */
/* SUBROUTINE SP2RLV(QHKL,QT,QM,QS,IER) */
/* SUBROUTINE INVS(S,SINV,IER) */
/* SUBROUTINE ERRESO(MODULE,IER) */
/* ------------------------------------------------------------------ */
/* File [MAD.SRC]T_RLP.FOR */
return 0;
} /* t_rlp__ */
/* Subroutine */ int setrlp_(real *sam, integer *ier)
{
/* System generated locals */
doublereal d__1;
/* Builtin functions */
double cos(doublereal), sin(doublereal), sqrt(doublereal), atan(
doublereal);
/* Local variables */
static doublereal alfa[3], cosa[3], cosb[3];
static integer imod;
static doublereal sina[3], sinb[3], aspv[6] /* was [3][2] */;
extern /* Subroutine */ int invs_(doublereal *, doublereal *, integer *);
static doublereal a[3], b[3], c__[3], bb[9] /* was [3][3] */, cc;
static integer id, ie, jd, je, jf, kg, lf, lh, md, me, ne;
static doublereal zp, vv[9] /* was [3][3] */;
extern /* Subroutine */ int erreso_(integer *, integer *);
static doublereal rlb[6] /* was [3][2] */;
/* ============================ */
/* SETRLP: Computation of matrix S which transforms (QH,QK,QL) to */
/* vector (Q1,Q2) in scattering plane (defined by vectors A1,A2) */
/* and SINV matrix for the inverse transformation */
/* INPUT SAM SAMPLE CHARACTERISTICS */
/* SAM(1)=AS LATTICE PARAMETERS */
/* SAM(2)=BS ------------------ */
/* SAM(3)=CS ------------------ */
/* SAM(4)=AA LATTICE ANGLES */
/* SAM(5)=BB -------------- */
/* SAM(6)=CC -------------- */
/* SAM(7)=AX VECTOR A IN SCATTERING PLANE */
/* SAM(8)=AY ---------------------------- */
/* SAM(9)=AZ ---------------------------- */
/* SAM(10)=BX VECTOR B IN SCATTERING PLANE */
/* SAM(11)=BY ---------------------------- */
/* SAM(12)=BZ ---------------------------- */
/* OUTPUT IER ERROR RETURN TO BE TREATED BY ERRESO */
/* IER=1 ERROR ON LATTICE PARAMETERS */
/* IER=2 ERROR ON LATTICE ANGLES */
/* IER=3 ERROR ON VECTORS A1, A2 */
/* ------------------------------------------------------------------ */
/* Part of [MAD.SRC]T_RLP.FOR */
/* ------------------------------------------------------------------ */
/* Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* DO NOT EXPORT THE FOLLOWING COMMON ! */
/* IT IS JUST FOR PERMANENT STORAGE USE */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* SOME TESTS AND INIT OF CALCUALTION */
/* Parameter adjustments */
--sam;
/* Function Body */
*ier = 0;
imod = 1;
zp = 6.2831853071795862;
osolem_1.iok = 0;
for (id = 1; id <= 3; ++id) {
a[id - 1] = sam[id];
alfa[id - 1] = sam[id + 3];
aspv[id - 1] = sam[id + 6];
aspv[id + 2] = sam[id + 9];
}
for (id = 1; id <= 3; ++id) {
*ier = 1;
if ((d__1 = a[id - 1], abs(d__1)) <= 1e-8) {
goto L999;
}
*ier = 0;
}
for (id = 1; id <= 3; ++id) {
a[id - 1] /= zp;
alfa[id - 1] /= 57.29577951308232087679815481410517;
cosa[id - 1] = cos(alfa[id - 1]);
sina[id - 1] = sin(alfa[id - 1]);
}
cc = cosa[0] * cosa[0] + cosa[1] * cosa[1] + cosa[2] * cosa[2];
cc = cosa[0] * 2. * cosa[1] * cosa[2] + 1. - cc;
*ier = 2;
if (cc <= .1) {
goto L999;
}
*ier = 0;
cc = sqrt(cc);
je = 2;
kg = 3;
for (id = 1; id <= 3; ++id) {
b[id - 1] = sina[id - 1] / (a[id - 1] * cc);
cosb[id - 1] = (cosa[je - 1] * cosa[kg - 1] - cosa[id - 1]) / (sina[
je - 1] * sina[kg - 1]);
sinb[id - 1] = sqrt(1. - cosb[id - 1] * cosb[id - 1]);
rlb[id + 2] = (d__1 = atan(sinb[id - 1] / cosb[id - 1]), abs(d__1)) *
57.29577951308232087679815481410517;
je = kg;
kg = id;
}
bb[0] = b[0];
bb[1] = 0.;
bb[2] = 0.;
bb[3] = b[1] * cosb[2];
bb[4] = b[1] * sinb[2];
bb[5] = 0.;
bb[6] = b[2] * cosb[1];
bb[7] = -b[2] * sinb[1] * cosa[0];
bb[8] = 1. / a[2];
for (id = 1; id <= 3; ++id) {
rlb[id - 1] = 0.;
for (je = 1; je <= 3; ++je) {
/* Computing 2nd power */
d__1 = bb[je + id * 3 - 4];
rlb[id - 1] += d__1 * d__1;
}
*ier = 1;
if ((d__1 = rlb[id - 1], abs(d__1)) <= 1e-8) {
goto L999;
}
*ier = 0;
rlb[id - 1] = sqrt(rlb[id - 1]);
}
/* ----------------------------------------------------------------------- */
/* GENERATION OF S ORIENTATION MATRIX REC. LATTICE TO SCATTERING PLANE */
for (kg = 1; kg <= 2; ++kg) {
for (ie = 1; ie <= 3; ++ie) {
vv[kg + ie * 3 - 4] = 0.;
for (jf = 1; jf <= 3; ++jf) {
vv[kg + ie * 3 - 4] += bb[ie + jf * 3 - 4] * aspv[jf + kg * 3
- 4];
}
}
}
for (md = 3; md >= 2; --md) {
for (ne = 1; ne <= 3; ++ne) {
id = md % 3 + 1;
je = (md + 1) % 3 + 1;
kg = ne % 3 + 1;
lh = (ne + 1) % 3 + 1;
vv[md + ne * 3 - 4] = vv[id + kg * 3 - 4] * vv[je + lh * 3 - 4] -
vv[id + lh * 3 - 4] * vv[je + kg * 3 - 4];
}
}
for (id = 1; id <= 3; ++id) {
c__[id - 1] = 0.;
for (je = 1; je <= 3; ++je) {
/* Computing 2nd power */
d__1 = vv[id + je * 3 - 4];
c__[id - 1] += d__1 * d__1;
}
*ier = 3;
if ((d__1 = c__[id - 1], abs(d__1)) <= 1e-6) {
goto L999;
}
*ier = 0;
c__[id - 1] = sqrt(c__[id - 1]);
}
for (id = 1; id <= 3; ++id) {
for (je = 1; je <= 3; ++je) {
vv[je + id * 3 - 4] /= c__[je - 1];
}
}
for (kg = 1; kg <= 3; ++kg) {
for (me = 1; me <= 3; ++me) {
osolem_1.s[kg + (me << 2) - 5] = 0.;
for (lf = 1; lf <= 3; ++lf) {
osolem_1.s[kg + (me << 2) - 5] += vv[kg + lf * 3 - 4] * bb[lf
+ me * 3 - 4];
}
}
}
osolem_1.s[15] = 1.;
for (jd = 1; jd <= 3; ++jd) {
osolem_1.s[(jd << 2) - 1] = 0.;
osolem_1.s[jd + 11] = 0.;
}
/* ----------------------------------------------------------------------- */
/* INVERT TRANSFORMATION MATRIX S AND PU RESULT IN SINV */
*ier = 3;
invs_(osolem_1.s, osolem_1.sinv, ier);
*ier = 0;
if (*ier != 0) {
goto L999;
}
osolem_1.iok = 123;
/* --------------------------------------------------------------------------- */
/* SORTIE */
L999:
if (*ier != 0) {
erreso_(&imod, ier);
}
return 0;
} /* setrlp_ */
/* Subroutine */ int rl2spv_(doublereal *qhkl, doublereal *qt, doublereal *qm,
doublereal *qs, integer *ier)
{
/* System generated locals */
doublereal d__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static integer id, je;
/* ========================================= */
/* ------------------------------------------------------------------ */
/* INPUT QHKL QHKL -> QT */
/* A Q-VECTOR TO BE TRANSFORM FROM RECIP LATTICE TO SCATTERING PLANE */
/* CHECK THAT Q-VECTOR IS IN THE PLANE */
/* INPUT Q-VECTOR QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE */
/* OUTPUT Q-VECTOR QT(3) Q-VECTOR IN SCATTERING PLANE */
/* OUTPUT QM AND QS QMODULUS AND ITS SQUARE ( TO BE VERIFIED ) */
/* OUTPUT ERROR IER */
/* IER=1 MATRIX S NOT OK */
/* IER=2 Q NOT IN SCATTERING PLANE */
/* IER=3 Q MODULUS TOO SMALL */
/* ------------------------------------------------------------------ */
/* Part of [MAD.SRC]T_RLP.FOR */
/* ------------------------------------------------------------------ */
/* Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* DO NOT EXPORT THE FOLLWING COOMON ! */
/* IT IS JUST FOR PERMANENT STORAGE USE */
/* --------------------------------------------------------------------------- */
/* --------------------------------------------------------------------------- */
/* INIT AND TEST IF TRANSFO MATRICES ARE OK */
/* Parameter adjustments */
--qt;
--qhkl;
/* Function Body */
*ier = 1;
if (osolem_1.iok != 123) {
goto L999;
}
*ier = 0;
/* ----------------------------------------------------------------------- */
for (id = 1; id <= 3; ++id) {
qt[id] = 0.;
for (je = 1; je <= 3; ++je) {
qt[id] += qhkl[je] * osolem_1.s[id + (je << 2) - 5];
}
}
*ier = 2;
if (abs(qt[3]) > 1e-4) {
goto L999;
}
*ier = 0;
*qs = 0.;
for (id = 1; id <= 3; ++id) {
/* Computing 2nd power */
d__1 = qt[id];
*qs += d__1 * d__1;
}
if (*qs < 1e-8) {
*ier = 3;
} else {
*qm = sqrt(*qs);
}
/* --------------------------------------------------------------------------- */
L999:
return 0;
} /* rl2spv_ */
/* Subroutine */ int sp2rlv_(doublereal *qhkl, doublereal *qt, doublereal *qm,
doublereal *qs, integer *ier)
{
/* System generated locals */
doublereal d__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static integer id, je;
/* ========================================= */
/* ------------------------------------------------------------------ */
/* INPUT QT QHKL <- QT */
/* A Q-VECTOR TO BE TRANSFORM FROM SCATTERING PLANE TO RECIP LATTICE */
/* CHECK THAT Q, D & G VECTORS ARE IN THE SCATTERING PLANE */
/* INPUT Q-VECTOR QT(3) Q-VECTOR IN SCATTERING PLANE */
/* OUTPUT Q-VECTOR QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE */
/* OUTPUT QM AND QS QMODULUS AND ITS SQUARE ( TO BE VERIFIED ) */
/* OUTPUT ERROR IER */
/* IER=1 MATRIX S NOT OK */
/* IER=2 Q NOT IN SCATTERING PLANE */
/* IER=3 Q MODULUS TOO SMALL */
/* ------------------------------------------------------------------ */
/* Part of [MAD.SRC]T_RLP.FOR */
/* ------------------------------------------------------------------ */
/* Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* DO NOT EXPORT THE FOLLWING COOMON ! */
/* IT IS JUST FOR PERMANENT STORAGE USE */
/* --------------------------------------------------------------------------- */
/* --------------------------------------------------------------------------- */
/* INIT AND TEST IF TRANSFO MATRICES ARE OK */
/* Parameter adjustments */
--qt;
--qhkl;
/* Function Body */
*ier = 1;
if (osolem_1.iok != 123) {
goto L999;
}
*ier = 2;
if (abs(qt[3]) > 1e-4) {
goto L999;
}
*ier = 0;
/* ----------------------------------------------------------------------- */
*qs = 0.;
for (id = 1; id <= 3; ++id) {
/* Computing 2nd power */
d__1 = qt[id];
*qs += d__1 * d__1;
}
if (*qs < 1e-8) {
*ier = 3;
} else {
*qm = sqrt(*qs);
}
/* ----------------------------------------------------------------------- */
for (id = 1; id <= 3; ++id) {
qhkl[id] = 0.;
for (je = 1; je <= 3; ++je) {
qhkl[id] += osolem_1.sinv[id + (je << 2) - 5] * qt[je];
}
}
/* --------------------------------------------------------------------------- */
L999:
return 0;
} /* sp2rlv_ */
/* Subroutine */ int invs_(doublereal *s, doublereal *sinv, integer *ier)
{
/* Initialized data */
static integer m[3] = { 2,3,1 };
static integer n[3] = { 3,1,2 };
static integer id, je, mi, mj, ni, nj;
static doublereal det;
/* ============================== */
/* ------------------------------------------------------------------ */
/* ROUTINE TO INVERT MATRIX S, GENERATED BY SETRLP, WHICH TRANSFORMS */
/* (QH,QK,QL) TO (Q1,Q2) IN THE SCATTERING PLANE */
/* INPUT MATRIX DOUBLE PRECISION S(4,4) */
/* OUTPUT MATRIX DOUBLE PRECISION SINV(4,4) */
/* OUTPUT ERROR IER */
/* IER=1 DETERMINANT OF MATRIX S TOO SMALL */
/* ------------------------------------------------------------------ */
/* Part of [MAD.SRC]T_RLP.FOR */
/* ------------------------------------------------------------------ */
/* Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
/* Parameter adjustments */
sinv -= 5;
s -= 5;
/* Function Body */
*ier = 0;
for (id = 1; id <= 4; ++id) {
for (je = 1; je <= 4; ++je) {
sinv[id + (je << 2)] = 0.;
}
}
det = 0.;
for (id = 1; id <= 3; ++id) {
for (je = 1; je <= 3; ++je) {
mi = m[id - 1];
mj = m[je - 1];
ni = n[id - 1];
nj = n[je - 1];
sinv[je + (id << 2)] = s[mi + (mj << 2)] * s[ni + (nj << 2)] - s[
ni + (mj << 2)] * s[mi + (nj << 2)];
}
det += s[id + 4] * sinv[(id << 2) + 1];
}
if (abs(det) < 1e-6) {
*ier = 1;
} else {
for (id = 1; id <= 3; ++id) {
for (je = 1; je <= 3; ++je) {
sinv[id + (je << 2)] /= det;
}
}
}
sinv[20] = 1.;
return 0;
} /* invs_ */
/* Subroutine */ int erreso_(integer *module, integer *ier)
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer lier, lmodule;
/* ============================= */
/* ------------------------------------------------------------------ */
/* SUBROUTINE TO TREAT ERRORS FROM RESOLUTION CALCULATIONS */
/* MODULE = 1 -> SETRLP */
/* MODULE = 2 -> RL2SPV */
/* MODULE = 3 -> EX_CASE */
/* ------------------------------------------------------------------ */
/* Part of [MAD.SRC]T_RLP.FOR */
/* include 'iolsddef.inc' */
/* ------------------------------------------------------------------ */
/* Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* --------------------------------------------------------------------------- */
/* Computing MIN */
i__1 = max(*ier,1);
lier = min(i__1,4);
/* Computing MIN */
i__1 = max(*module,1);
lmodule = min(i__1,3);
/* WRITE(iolun,501) MESER(LIER,LMODULE) */
/* L501: */
return 0;
} /* erreso_ */

463
t_rlp.f Executable file
View File

@ -0,0 +1,463 @@
SUBROUTINE T_RLP ! File [MAD.SRC]T_RLP.FOR
c ================
c
cdec$ Ident 'V01A'
c------------------------------------------------------------------
c Updates:
c V01A 7-May-1996 DM. Put error output to IOLUN, use IMPLICIT NONE and
c get the code indented so that it is readable!
c------------------------------------------------------------------
c Routines to deal with the reciprocical lattice PB
c------------------------------------------------------------------
c Entry points in this file:
c
c SETRLP : CALCULATION OF S AND INVS , ORIENTATION MATRIX
C RL2SPV : TRANSFO FROM RECIP LAT TO SCAT PLANE
C SP2RLV : TRANSFO FROM SCAT PLANE TO RECIP LAT
C INVS : INVERT MATRIX S, GENERATED BY SETRLP.
C ERRESO : DEAL ITH ERROR MESSAGES FOR ALL MODULES
C
C SUBROUTINE SETRLP(SAM,IER)
C SUBROUTINE RL2SPV(QHKL,QT,QM,QS,IER)
C SUBROUTINE SP2RLV(QHKL,QT,QM,QS,IER)
C SUBROUTINE INVS(S,SINV,IER)
C SUBROUTINE ERRESO(MODULE,IER)
c------------------------------------------------------------------
implicit none
end
c
SUBROUTINE SETRLP (SAM, IER) ! Part of [MAD.SRC]T_RLP.FOR
c ============================
c
C SETRLP: Computation of matrix S which transforms (QH,QK,QL) to
C vector (Q1,Q2) in scattering plane (defined by vectors A1,A2)
C and SINV matrix for the inverse transformation
C
C INPUT SAM SAMPLE CHARACTERISTICS
C SAM(1)=AS LATTICE PARAMETERS
C SAM(2)=BS ------------------
C SAM(3)=CS ------------------
C SAM(4)=AA LATTICE ANGLES
C SAM(5)=BB --------------
C SAM(6)=CC --------------
C SAM(7)=AX VECTOR A IN SCATTERING PLANE
C SAM(8)=AY ----------------------------
C SAM(9)=AZ ----------------------------
C SAM(10)=BX VECTOR B IN SCATTERING PLANE
C SAM(11)=BY ----------------------------
C SAM(12)=BZ ----------------------------
C OUTPUT IER ERROR RETURN TO BE TREATED BY ERRESO
C IER=1 ERROR ON LATTICE PARAMETERS
C IER=2 ERROR ON LATTICE ANGLES
C IER=3 ERROR ON VECTORS A1, A2
c------------------------------------------------------------------
IMPLICIT NONE
c
DOUBLE PRECISION PI
PARAMETER (PI=3.14159265358979323846264338327950D0)
DOUBLE PRECISION RD
PARAMETER (RD=57.29577951308232087679815481410517D0)
DOUBLE PRECISION EPS
PARAMETER (EPS=1.D-1)
DOUBLE PRECISION EPS6
PARAMETER (EPS6=1.D-6)
DOUBLE PRECISION EPS8
PARAMETER (EPS8=1.D-8)
c------------------------------------------------------------------
C Define the dummy arguments
real*4 sam(64)
integer*4 ier
c------------------------------------------------------------------
C DO NOT EXPORT THE FOLLOWING COMMON !
C IT IS JUST FOR PERMANENT STORAGE USE
C
DOUBLE PRECISION S, SINV
integer*4 iok
COMMON /OSOLEM/ S(4,4), SINV(4,4), iok
C-----------------------------------------------------------------------
double precision a(3), aspv(3,2), rlb(3,2)
double precision alfa(3), sina(3), sinb(3), cosa(3), cosb(3)
double precision b(3), c(3)
double precision vv(3,3), bb(3,3)
C
double precision zp, cc
integer*4 imod
integer*4 id, ie
integer*4 jd, je, jf
integer*4 kg
integer*4 lh, lf
integer*4 md, me
integer*4 ne
C-----------------------------------------------------------------------
C SOME TESTS AND INIT OF CALCUALTION
C
ier = 0
IMOD = 1
ZP = 2.D0*PI
IOK = 0
DO ID = 1,3
A(ID) = SAM(ID)
ALFA(ID) = SAM(ID+3)
ASPV(ID,1) = SAM(6+ID)
ASPV(ID,2) = SAM(9+ID)
ENDDO
C
DO ID = 1,3
IER = 1
IF(ABS(A(ID)).LE.EPS8) GOTO 999
IER = 0
ENDDO
DO ID = 1,3
A(ID) = A(ID)/ZP
ALFA(ID) = ALFA(ID)/RD
COSA(ID) = COS(ALFA(ID))
SINA(ID) = SIN(ALFA(ID))
ENDDO
CC = COSA(1)*COSA(1)+COSA(2)*COSA(2)+COSA(3)*COSA(3)
CC = (1.D0+2.D0*COSA(1)*COSA(2)*COSA(3)-CC)
IER = 2
IF(CC.LE.EPS) GOTO 999
IER = 0
CC = SQRT(CC)
JE = 2
KG = 3
DO ID = 1,3
B(ID) = SINA(ID)/(A(ID)*CC)
COSB(ID) = (COSA(JE)*COSA(KG)-COSA(ID))/(SINA(JE)*SINA(KG))
SINB(ID) = SQRT(1.D0-COSB(ID)*COSB(ID))
RLB(ID,2) = ABS(ATAN(SINB(ID)/COSB(ID)))*RD
JE = KG
KG = ID
ENDDO
BB(1,1) = B(1)
BB(2,1) = 0.D0
BB(3,1) = 0.D0
BB(1,2) = B(2)*COSB(3)
BB(2,2) = B(2)*SINB(3)
BB(3,2) = 0.D0
BB(1,3) = B(3)*COSB(2)
BB(2,3) = -B(3)*SINB(2)*COSA(1)
BB(3,3) = 1.D0/A(3)
C
DO ID = 1,3
RLB(ID,1) = 0.D0
DO JE = 1,3
RLB(ID,1) = RLB(ID,1)+BB(JE,ID)**2
ENDDO
IER = 1
IF (ABS(RLB(ID,1)).LE.EPS8) GOTO 999
IER = 0
RLB(ID,1) = SQRT(RLB(ID,1))
ENDDO
C-----------------------------------------------------------------------
C GENERATION OF S ORIENTATION MATRIX REC. LATTICE TO SCATTERING PLANE
C
DO KG = 1,2
DO IE = 1,3
VV(KG,IE) = 0.D0
DO JF = 1,3
VV(KG,IE) = VV(KG,IE)+BB(IE,JF)*ASPV(JF,KG)
ENDDO
ENDDO
ENDDO
DO MD = 3,2,-1
DO NE = 1,3
ID = MOD(MD,3)+1
JE = MOD(MD+1,3)+1
KG = MOD(NE,3)+1
LH = MOD(NE+1,3)+1
VV(MD,NE) = VV(ID,KG)*VV(JE,LH)-VV(ID,LH)*VV(JE,KG)
ENDDO
ENDDO
C
DO ID = 1,3
C(ID) = 0.D0
DO JE = 1,3
C(ID) = C(ID)+VV(ID,JE)**2
ENDDO
IER = 3
IF (ABS(C(ID)).LE.EPS6) GOTO 999
IER = 0
C(ID) = SQRT(C(ID))
ENDDO
C
DO ID = 1,3
DO JE = 1,3
VV(JE,ID) = VV(JE,ID)/C(JE)
ENDDO
ENDDO
DO KG = 1,3
DO ME = 1,3
S(KG,ME) = 0.D0
DO LF = 1,3
S(KG,ME) = S(KG,ME)+VV(KG,LF)*BB(LF,ME)
ENDDO
ENDDO
ENDDO
S(4,4) = 1.D0
DO JD = 1,3
S(4,JD) = 0.D0
S(JD,4) = 0.D0
ENDDO
C-----------------------------------------------------------------------
C INVERT TRANSFORMATION MATRIX S AND PU RESULT IN SINV
C
IER = 3
CALL INVS(S,SINV,IER)
IER = 0
IF (IER.NE.0)GOTO 999
IOK = 123
C---------------------------------------------------------------------------
C SORTIE
C
999 CONTINUE
IF (IER.NE.0) CALL ERRESO(IMOD,IER)
RETURN
END
c
SUBROUTINE RL2SPV (QHKL, QT, QM, QS, IER) ! Part of [MAD.SRC]T_RLP.FOR
c =========================================
c
c------------------------------------------------------------------
C INPUT QHKL QHKL -> QT
C A Q-VECTOR TO BE TRANSFORM FROM RECIP LATTICE TO SCATTERING PLANE
C CHECK THAT Q-VECTOR IS IN THE PLANE
C
C INPUT Q-VECTOR QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE
C
C OUTPUT Q-VECTOR QT(3) Q-VECTOR IN SCATTERING PLANE
C OUTPUT QM AND QS QMODULUS AND ITS SQUARE ( TO BE VERIFIED )
C OUTPUT ERROR IER
C IER=1 MATRIX S NOT OK
C IER=2 Q NOT IN SCATTERING PLANE
C IER=3 Q MODULUS TOO SMALL
c------------------------------------------------------------------
IMPLICIT NONE
c
DOUBLE PRECISION DEPS4
PARAMETER (DEPS4=1.D-4)
DOUBLE PRECISION DEPS8
PARAMETER (DEPS8=1.D-8)
c------------------------------------------------------------------
c Define the dummy arguments
DOUBLE PRECISION QHKL(3)
DOUBLE PRECISION QT(3)
DOUBLE PRECISION QM
DOUBLE PRECISION QS
integer*4 IER
c------------------------------------------------------------------
C DO NOT EXPORT THE FOLLWING COOMON !
C IT IS JUST FOR PERMANENT STORAGE USE
C
DOUBLE PRECISION S, SINV
integer*4 iok
COMMON /OSOLEM/ S(4,4), SINV(4,4), iok
C---------------------------------------------------------------------------
integer*4 id, je
C---------------------------------------------------------------------------
C INIT AND TEST IF TRANSFO MATRICES ARE OK
IER = 1
IF (IOK.NE.123) GOTO 999
IER = 0
C-----------------------------------------------------------------------
C
DO ID = 1,3
QT(ID) = 0.D0
DO JE = 1,3
QT(ID) = QT(ID) + QHKL(JE)*S(ID,JE)
ENDDO
ENDDO
IER = 2
IF(ABS(QT(3)).GT.DEPS4) GOTO 999
IER = 0
QS = 0.D0
DO ID = 1,3
QS = QS+QT(ID)**2
ENDDO
IF(QS.LT.DEPS8) THEN
IER = 3
ELSE
QM = SQRT(QS)
ENDIF
C---------------------------------------------------------------------------
C
999 CONTINUE
RETURN
END
c
SUBROUTINE SP2RLV (QHKL, QT, QM, QS, IER) ! Part of [MAD.SRC]T_RLP.FOR
c =========================================
c
c------------------------------------------------------------------
C INPUT QT QHKL <- QT
C A Q-VECTOR TO BE TRANSFORM FROM SCATTERING PLANE TO RECIP LATTICE
C CHECK THAT Q, D & G VECTORS ARE IN THE SCATTERING PLANE
C
C INPUT Q-VECTOR QT(3) Q-VECTOR IN SCATTERING PLANE
C
C OUTPUT Q-VECTOR QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE
C OUTPUT QM AND QS QMODULUS AND ITS SQUARE ( TO BE VERIFIED )
C OUTPUT ERROR IER
C IER=1 MATRIX S NOT OK
C IER=2 Q NOT IN SCATTERING PLANE
C IER=3 Q MODULUS TOO SMALL
c------------------------------------------------------------------
IMPLICIT NONE
C
DOUBLE PRECISION EPS4
PARAMETER (EPS4=1.D-4)
DOUBLE PRECISION EPS8
PARAMETER (EPS8=1.D-8)
c------------------------------------------------------------------
c Define the dummy arguments
DOUBLE PRECISION QHKL(3)
DOUBLE PRECISION QT(3)
DOUBLE PRECISION QM
DOUBLE PRECISION QS
integer*4 IER
c------------------------------------------------------------------
C DO NOT EXPORT THE FOLLWING COOMON !
C IT IS JUST FOR PERMANENT STORAGE USE
C
DOUBLE PRECISION S, SINV
integer*4 iok
COMMON /OSOLEM/ S(4,4), SINV(4,4), iok
C---------------------------------------------------------------------------
integer*4 id, je
C---------------------------------------------------------------------------
C INIT AND TEST IF TRANSFO MATRICES ARE OK
IER = 1
IF (IOK.NE.123) GOTO 999
IER = 2
IF(ABS(QT(3)).GT.EPS4) GOTO 999
IER = 0
C-----------------------------------------------------------------------
QS = 0.D0
DO ID = 1,3
QS = QS+QT(ID)**2
ENDDO
IF(QS.LT.EPS8) THEN
IER = 3
ELSE
QM = SQRT(QS)
ENDIF
C-----------------------------------------------------------------------
C
DO ID = 1,3
QHKL(ID) = 0.D0
DO JE = 1,3
QHKL(ID) = QHKL(ID)+SINV(ID,JE)*QT(JE)
ENDDO
ENDDO
C---------------------------------------------------------------------------
C
999 CONTINUE
RETURN
END
c
SUBROUTINE INVS (S, SINV, IER) ! Part of [MAD.SRC]T_RLP.FOR
c ==============================
c
c------------------------------------------------------------------
C ROUTINE TO INVERT MATRIX S, GENERATED BY SETRLP, WHICH TRANSFORMS
C (QH,QK,QL) TO (Q1,Q2) IN THE SCATTERING PLANE
C INPUT MATRIX DOUBLE PRECISION S(4,4)
C OUTPUT MATRIX DOUBLE PRECISION SINV(4,4)
C OUTPUT ERROR IER
C IER=1 DETERMINANT OF MATRIX S TOO SMALL
c------------------------------------------------------------------
IMPLICIT NONE
c
DOUBLE PRECISION EPS6
PARAMETER (EPS6=1.D-6)
c------------------------------------------------------------------
c Define the dummy arguments
DOUBLE PRECISION S(4,4)
DOUBLE PRECISION SINV(4,4)
integer*4 IER
c------------------------------------------------------------------
integer*4 m(3)
integer*4 n(3)
c
integer*4 id, je
integer*4 mi, mj, ni, nj
double precision det
c------------------------------------------------------------------
DATA M/2,3,1/
DATA N/3,1,2/
IER = 0
DO ID = 1,4
DO JE = 1,4
SINV(ID,JE) = 0.D0
ENDDO
ENDDO
DET = 0.D0
DO ID = 1,3
DO JE = 1,3
MI = M(ID)
MJ = M(JE)
NI = N(ID)
NJ = N(JE)
SINV(JE,ID) = S(MI,MJ)*S(NI,NJ)-S(NI,MJ)*S(MI,NJ)
ENDDO
DET = DET+S(ID,1)*SINV(1,ID)
ENDDO
IF(ABS(DET).LT.EPS6) THEN
IER = 1
ELSE
DO ID = 1,3
DO JE = 1,3
SINV(ID,JE) = SINV(ID,JE)/DET
ENDDO
ENDDO
ENDIF
SINV(4,4) = 1.D0
RETURN
END
c
SUBROUTINE ERRESO(MODULE,IER) ! Part of [MAD.SRC]T_RLP.FOR
c =============================
c
c------------------------------------------------------------------
C SUBROUTINE TO TREAT ERRORS FROM RESOLUTION CALCULATIONS
C MODULE = 1 -> SETRLP
C MODULE = 2 -> RL2SPV
C MODULE = 3 -> EX_CASE
c------------------------------------------------------------------
IMPLICIT NONE
c
integer*4 MXER
PARAMETER (MXER = 4)
integer*4 MXMOD
PARAMETER (MXMOD = 3)
c
C include 'iolsddef.inc'
c------------------------------------------------------------------
c Define the dummy arguments
integer*4 module
integer*4 ier
c------------------------------------------------------------------
integer*4 lmodule, lier
CHARACTER*64 MESER(MXER,MXMOD)
DATA MESER/
+ ' Check lattice spacings (AS,BS,CS)',
+ ' Check cell angles (AA,BB,CC)',
+ ' Check scattering plane (AX....BZ)',
+ ' Check lattice and scattering plane',
C
+ ' Check Lattice and Scattering Plane',
+ ' Q not in scattering plane',
+ ' Q modulus too small',
+ ' KI,KF,Q triangle cannot close',
C
+ ' Error in KI or KF. Check D-spacings & units',
+ ' KI or KF cannot be obtained',
+ ' KI or KF too small',
+ ' KI,KF,Q triangle will not close'/
C---------------------------------------------------------------------------
LIER = MIN(MAX(IER,1),MXER)
LMODULE = MIN(MAX(MODULE,1),MXMOD)
C WRITE(iolun,501) MESER(LIER,LMODULE)
501 FORMAT(A)
RETURN
END

495
t_update.c Normal file
View File

@ -0,0 +1,495 @@
/* 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__ */

409
t_update.f Executable file
View File

@ -0,0 +1,409 @@
SUBROUTINE T_UPDATE ( ! File [MAD.SRC]T_UPDATE.FOR
c ===================
+ 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)
c
cdec$ Ident 'V01C'
c------------------------------------------------------------------
c Note:
c IF1,If2 changed to Int*4. They were Real*4 in ILL version.
c------------------------------------------------------------------
c Updates:
c V01A 7-May-1996 DM. Put error output to IOLUN, use IMPLICIT NONE and
c get the code indented so that it is readable!
c V01C 12-Oct-1998 DM. Put in Mag Field calculation for SINQ Helmolz coils.
C-----------------------------------------------------------------------
c Entry points in this file:
C T_UPDATE : USE THE MOTOR ANGLES TO CALCULATE VALUES OF
c EI,KI,EF,KF,QH,QK,QL,EN AND TO DETERMINE
c WHETHER FLIPPERS F1,F2 ARE 'ON' OR 'OFF'.
c EX_UP : CALCULATE EX AND AKX FORM AX2 VALUES.
c SAM_UP : CALCULATE THINGS IN RECIPROCICAL LATTICE.
c HELM_UP : CALCULATE HELMHOLTZ FIELDS.
c FLIP_UP : CHECK IF FLIPPERS ON OR OFF.
C-----------------------------------------------------------------------
c T_UPDATE:
C ROUTINE TO USE THE MOTOR ANGLES TO
C CALCULATE VALUES OF EI,KI,EF,KF,QH,QK,QL,EN
C USE CURRENT-SUPPLY VALUES TO COMPUTE HELMHOLTZ FIELDS HX,HY,HZ
C AND TO DETERMINE WHETHER FLIPPERS F1,F2 ARE 'ON' OR 'OFF'.
C FLIPPERS ARE 'ON' ONLY IF CURRENTS ARE THOSE GIVEN BY
C F1V,F1H,F2V,F2H
C-----------------------------------------------------------------------
IMPLICIT NONE
ccc IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C include 'iolsddef.inc'
C-----------------------------------------------------------------------
c Define the dummy arguments
c Input:
real*4 p_a(6) ! Positions of angles A1-A6.
real*4 p_ih(8) ! Position of currents for flippers and ..
! .. helmotz (8 currents).
real*4 c_ih(4) ! Conversion factors for Helmotz (4 currents).
! Configuration of the machine ...
logical*4 lpa ! True if machine in polarization mode.
real*4 dm ! Monochr. d-spacing
real*4 da ! Analyser d-spacing
integer*4 isa ! Scattering sense at analyser (probably) ..
! .. +ve to the left.
real*4 helm ! Angle between axis of Helmolz pair one and ..
! .. ki (probably).
real*4 f1h, f1v ! Flipper 1 hor and vert currents (hor ..
! .. current is possibly ki*f1h ??)
real*4 f2h, f2v ! Flipper 2 hor and vert currents (hor ..
! .. current is possibly kf*f2h ??)
real*4 f ! Energy unit
c Output:
real*4 ei ! Incident neutron energy
real*4 aki ! Incident neutron wave vector
real*4 ef ! Final neutron energy
real*4 akf ! Final neutron wave vector
real*4 qhkl(3) ! Components of q in reciprocal space
real*4 en ! Energy transfer
real*4 hx, hy, hz ! Components of Helmolz field on sample
integer*4 if1, if2 ! Status of Flippers 1 and 2
real*4 qm ! Length of q
integer*4 ier ! Error status
C-----------------------------------------------------------------------
double precision dbqhkl(3), dhx, dhy, dhz
double precision ddm, dda, df
double precision da2, da3, da4, da6
double precision dei, daki, def, dakf
double precision dqm, dqs, dphi
real*4 qs
c
integer*4 ieri, ieru, imod, id
C-----------------------------------------------------------------------
C SET UP
C
DDM=DM
DDA=DA
DF=F
DA2=P_A(2)
DA3=P_A(3)
DA4=P_A(4)
DA6=P_A(6)
C-----------------------------------------------------------------------
C
IERI=0
IERU=1 + 8
CALL EX_UP(DDM,DEI,DAKI,DA2,DF,IERI)
IF (IERI.EQ.0) THEN
EI=DEI
AKI=DAKI
IERU=0
ELSE
IMOD=3
CALL ERRESO(IMOD,IERI)
IERU = IERI + 8
ENDIF
IERI=0
IERU=1 + 8
CALL EX_UP(DDA,DEF,DAKF,DA6,DF,IERI)
IF (IERI.EQ.0) THEN
EF=DEF
AKF=DAKF
IERU=0
ELSE
IF (ISA.EQ.0) THEN
EF=DEI
AKF=DAKI
DEF=DEI
DAKF=DAKI
IERU=0
ELSE
IMOD=3
CALL ERRESO(IMOD,IERI)
IERU = 8 + IERI
ENDIF
ENDIF
IF (ISA.EQ.0) THEN
EF=DEI
AKF=DAKI
DEF=DEI
DAKF=DAKI
IERU=0
ENDIF
IF (IERU.EQ.0) EN=DEI-DEF
IERI=0
IERU=1 + 4
CALL SAM_UP(DBQHKL,DQM,DQS,DPHI,DAKI,DAKF,DA3,DA4,IERI)
IF (IERI.EQ.0) THEN
DO ID=1,3
QHKL(ID)=DBQHKL(ID)
ENDDO
QM=DQM
QS=DQS
IERU=0
ELSE
IMOD=2
CALL ERRESO(IMOD,IERI)
IERU = IERI + 4
ENDIF
C
IERI=0
IF (LPA) THEN
IERU=1
CALL HELM_UP(DPHI,HELM,P_IH,C_IH,DHX,DHY,DHZ,IERI)
IF (IERI.NE.0) THEN
C WRITE(6,*) 'ERROR IN HELM_UP CHECK COEFF'
ELSE
HX=DHX
HY=DHY
HZ=DHZ
IERU=0
ENDIF
CALL FLIP_UP(IF1,IF2,P_IH,F1V,F1H,F2V,F2H,AKI,AKF,IERI)
ENDIF
C-----------------------------------------------------------------------
IER=IERU
RETURN
END
c
SUBROUTINE EX_UP(DX,EX,AKX,AX2,F,IER) ! Part of [MAD.SRC]T_UPDATE.FOR
c ================
c
C-----------------------------------------------------------------------
C CALCULATE EX AND AKX FORM AX2 VALUES
C-----------------------------------------------------------------------
implicit none
c
DOUBLE PRECISION PI
PARAMETER (PI=3.1415926535897932384626433832795E0)
DOUBLE PRECISION RD
PARAMETER (RD=57.29577951308232087679815481410517E0)
DOUBLE PRECISION EPS4
PARAMETER (EPS4=1.D-4)
C-----------------------------------------------------------------------
c Define the dummy arguments
c Input:
DOUBLE PRECISION DX ! D-SPACING OF CRISTAL
DOUBLE PRECISION AX2 ! TAKE OFF ANGLE
DOUBLE PRECISION F ! ENERGY UNIT
c Output:
DOUBLE PRECISION EX ! ENERGY
DOUBLE PRECISION AKX ! MOMENTUM
INTEGER*4 IER ! Error - IER=1 if DX OR AX TOO SMALL
C-----------------------------------------------------------------------
DOUBLE PRECISION ARG
C-----------------------------------------------------------------------
C!!!!!!!!!! This has to be fixed manually after conversion by f2c.
C!!!!!!!!!! The reason is a different definition of the abs function.
C!!!!!!!!!! MK, May 2001
ARG=ABS(DX*SIN(AX2/(2.D0*RD)))
IF(ARG.LE.EPS4) THEN
IER=1
ELSE
IER=0
AKX=PI/ARG
EX=F*AKX**2
ENDIF
C-----------------------------------------------------------------------
RETURN
END
c
SUBROUTINE SAM_UP ( ! Part of [MAD.SRC]T_UPDATE.FOR
c =================
+ QHKL, QM, QS, PHI,
+ AKI, AKF, A3, A4, IER)
C-----------------------------------------------------------------------
implicit none
c
double precision PI
PARAMETER (PI=3.1415926535897932384626433832795D0)
double precision RD
PARAMETER (RD=57.29577951308232087679815481410517D0)
double precision EPS3
PARAMETER (EPS3=1.D-3)
C-----------------------------------------------------------------------
c Define the dummy arguments
c Input:
double precision AKI ! Actual value of KI
double precision AKF ! Actual value of KF
double precision A3 ! Actual value of A3
double precision A4 ! Actual value of A4
c Output:
double precision QHKL(3) ! CALCULATED POSITIONS IN RECIP. LATTICE
double precision QM ! Q MODULUS
double precision QS ! Q MODULUS SQUARED
double precision PHI ! ANGLE BETWEEN KI AND Q
integer*4 IER ! Error status. IER=1 if QM TOO SMALL (PHI
! NOT CALCULATED)
C-----------------------------------------------------------------------
C LOCAL PARAMETERS
double precision QT(3)
double precision QPAR, QPERP
C-----------------------------------------------------------------------
C QPAR AND QPERP ARE COMPONENTS OF Q PARALLEL AND PERP TO KI
C TRANSFORM FROM SCATTERING PLANE INTO RECIPROCAL LATTICE
C
qpar = aki - akf * cos(a4/RD)
qperp = -akf * sin(a4/RD)
qt(1) = qpar * cos(a3/RD) + qperp * sin(a3/RD)
qt(2) = -qpar * sin(a3/RD) + qperp * cos(a3/RD)
qt(3) = 0.d0
call sp2rlv (qhkl, qt, qm, qs, ier)
ier = 3
if (qm .gt. EPS3) ier = 0
phi = 0.d0
if (abs(qperp) .gt. EPS3 .and. abs(qpar) .gt. EPS3) then
phi = atan2 (qperp, qpar)
elseif (abs(qpar) .lt. EPS3) then
if (a4 .gt. 0.0) then
phi = -0.5 * PI
else
phi = 0.5 * PI
endif
endif
c-----------------------------------------------------------------------
return
end
c
SUBROUTINE HELM_UP ( ! Part of [MAD.SRC]T_UPDATE.FOR
c ==================
+ PHI, HELM, P_IH, C_IH,
+ DHX, DHY, DHZ, IER)
c
C----------------------------------------------------------------
C HELMHOLTZ COILS UPDATE (ONLY IF LPA IS TRUE)
C----------------------------------------------------------------
implicit none
c
C include 'common_sinq.inc'
c
double precision PI
PARAMETER (PI = 3.1415926535897932384626433832795D0)
double precision RD
PARAMETER (RD = 57.29577951308232087679815481410517D0)
double precision EPS3
PARAMETER (EPS3 = 1.0D-3)
C-----------------------------------------------------------------------
c At ILL:
c There are 3 coils for Hx/Hy at 120 degrees to each other.
c
c There is a 4th coil for Hz.
c
c At SINQ:
c There is an Hx coil and an Hy coil (actually each is 4 coils powered
c in series). They are mounted on a ring (SRO). The value of HELM is
c the angle between the Hx coil and ki.
c
c There is a 3rd coil for Hz.
C-----------------------------------------------------------------------
c Define the dummy arguments
c Input:
double precision PHI ! Angle between KI and Q in scatt'g plane (rad)
real*4 HELM ! Angle between first coil and KI (degrees)
real*4 P_IH(8) ! Current values of 4 currents associated ..
! with Helmolz coils
real*4 C_IH(4) ! Conversion factor between Amperes and Gauss
c Output:
double precision DHX ! \
double precision DHY ! > The calculated fields
double precision DHZ ! /
integer*4 IER ! Error status. IER=1 if C_I HAS A ZERO COEF
C----------------------------------------------------------------
C LOCAL VARIABLES
integer*4 ic
double precision h(4), hpar, hperp, hmod, hdir, hdir2
LOGICAL AT_SINQ
C----------------------------------------------------------------
C H !FIELD OF 4 COILS AROUND SAMPLE
C HPAR !FIELD PAR COIL 1
C HPERP !FIELD PERP. TO COIL 1
C HDIR2 !ANGLE BETWEEN FIELD AND COIL 1
C HDIR !ANGLE BETWEEN FIELD AND Q
C----------------------------------------------------------------
ier = 0
AT_SINQ = .TRUE.
if (.not. at_sinq) then
do ic = 1, 4
h(ic) = 0.d0
if (abs(c_ih(ic)) .gt. EPS3) then
h(ic) = p_ih(ic+4) * c_ih(ic)
else
ier = 1
endif
enddo
hpar = h(1) - .5d0 * (h(2) + h(3))
hperp = .5d0 * sqrt(3.d0) * (h(2) - h(3))
dhz = h(4)
else
do ic = 1, 3
h(ic) = 0.d0
if (abs(c_ih(ic)) .gt. EPS3) then
h(ic) = p_ih(ic+4) * c_ih(ic)
else
ier = 1
endif
enddo
hpar = h(1)
hperp = h(2)
dhz = h(3)
endif
c
hmod = sqrt (hpar**2 + hperp**2)
if (abs(hpar) .gt. EPS3 .and. abs(hperp) .gt. EPS3) then
hdir2 = atan2 (abs(hperp), abs(hpar))
if (hpar .gt. 0 .and. hperp .lt. 0) then
hdir2 = -hdir2
elseif (hpar .lt. 0 .and. hperp .gt. 0) then
hdir2 = PI - hdir2
elseif (hpar .lt. 0 .and. hperp .lt. 0) then
hdir2 = hdir2 - PI
endif
elseif (abs(hpar) .gt. EPS3) then
if (hpar .ge. 0.0) hdir2 = 0.0
if (hpar .lt. 0.0) hdir2 = PI
elseif (abs(hperp) .gt. EPS3) then
if (hperp .ge. 0.0) hdir2 = 0.5 * PI
if (hperp .lt. 0.0) hdir2 = -0.5 * PI
else
hdir2 = 0.0
endif
hdir = hdir2 + (helm/RD) - phi
dhx = hmod * cos(hdir)
dhy = hmod * sin(hdir)
c
return
end
c
SUBROUTINE FLIP_UP ( ! Part of [MAD.SRC]T_UPDATE.FOR
c ==================
+ IF1, IF2,
+ P_IH, F1V, F1H, F2V, F2H,
+ AKI, AKF, IER)
C----------------------------------------------------------------
C FLIPPERS; ONLY IF LPA
C----------------------------------------------------------------
implicit none
C-----------------------------------------------------------------------
c Define the dummy arguments
c Input:
real*4 p_ih(8) ! Current values of 4 currents associated ..
! with flipper coils
real*4 f1v, f1h ! Flipper 1 vert and hor currents (hor ..
! .. current is possibly aki*f1h ??)
real*4 f2v, f2h ! Flipper 2 vert and hor currents (hor ..
! .. current is possibly akf*f2h ??)
real*4 aki ! Incident neutron wave vector
real*4 akf ! Final neutron wave vector
c Output:
integer*4 if1 ! Status of Flipper 1 (=1 if on and set OK)
integer*4 if2 ! Status of Flipper 2 (=1 if on and set OK)
integer*4 ier ! Error status. 0 if no error.
C----------------------------------------------------------------
IER = 0
IF1 = 0
IF2 = 0
IF ((ABS(P_IH(1) - F1V) .LT. 0.05) .AND.
+ (ABS(P_IH(2) - AKI*F1H) .LT. 0.05)) IF1 = 1
IF ((ABS(P_IH(3) - F2V) .LT. 0.05) .AND.
+ (ABS(P_IH(4) - AKF*F2H) .LT. 0.05)) IF2 = 1
C-----------------------------------------------------------------------
RETURN
END