From cd13637987d280230deccc219f654bd2f99e1387 Mon Sep 17 00:00:00 2001 From: cvs Date: Mon, 30 Jun 2003 11:51:39 +0000 Subject: [PATCH] - Updated Makefiles - Moved TAS code to psi - Updated programmers documentation --- hardsup/el734_utility.c | 9 + hardsup/makefile_alpha | 2 +- make_gen | 3 +- nextrics.c | 17 +- t_conv.c | 838 ++++++++++++++++++++++++++++++++++++++++ t_conv.f | 690 +++++++++++++++++++++++++++++++++ t_rlp.c | 516 +++++++++++++++++++++++++ t_rlp.f | 463 ++++++++++++++++++++++ t_update.c | 495 ++++++++++++++++++++++++ t_update.f | 409 ++++++++++++++++++++ 10 files changed, 3439 insertions(+), 3 deletions(-) create mode 100644 t_conv.c create mode 100755 t_conv.f create mode 100644 t_rlp.c create mode 100755 t_rlp.f create mode 100644 t_update.c create mode 100755 t_update.f diff --git a/hardsup/el734_utility.c b/hardsup/el734_utility.c index 53efe51..acdfbb7 100644 --- a/hardsup/el734_utility.c +++ b/hardsup/el734_utility.c @@ -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; diff --git a/hardsup/makefile_alpha b/hardsup/makefile_alpha index 7cd80b9..3c58898 100644 --- a/hardsup/makefile_alpha +++ b/hardsup/makefile_alpha @@ -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 diff --git a/make_gen b/make_gen index ff27176..116e5a5 100644 --- a/make_gen +++ b/make_gen @@ -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 diff --git a/nextrics.c b/nextrics.c index 6aa9bc2..bfa7db5 100644 --- a/nextrics.c +++ b/nextrics.c @@ -561,6 +561,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 */ @@ -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 diff --git a/t_conv.c b/t_conv.c new file mode 100644 index 0000000..6e6c580 --- /dev/null +++ b/t_conv.c @@ -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__ */ + diff --git a/t_conv.f b/t_conv.f new file mode 100755 index 0000000..ba39ef3 --- /dev/null +++ b/t_conv.f @@ -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 diff --git a/t_rlp.c b/t_rlp.c new file mode 100644 index 0000000..8d66baa --- /dev/null +++ b/t_rlp.c @@ -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_ */ + diff --git a/t_rlp.f b/t_rlp.f new file mode 100755 index 0000000..26b7b69 --- /dev/null +++ b/t_rlp.f @@ -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 diff --git a/t_update.c b/t_update.c new file mode 100644 index 0000000..99b9655 --- /dev/null +++ b/t_update.c @@ -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__ */ + diff --git a/t_update.f b/t_update.f new file mode 100755 index 0000000..f5c63d8 --- /dev/null +++ b/t_update.f @@ -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