691 lines
22 KiB
Fortran
Executable File
691 lines
22 KiB
Fortran
Executable File
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: (Vertical)
|
||
c CM1RX + CM2RX/SIN(abs(A1)/RD)
|
||
c
|
||
c For analyser: (Horizontal)
|
||
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
|