/* 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__ */