Files
sics/tacov.f

432 lines
16 KiB
Fortran

C-----------------------------------------------------------------------
C FILE T_CONV
C SUBROUTINE T_CONV(EI,AKI,EF,AKF,QHKL,EN,HX,HY,HZ,IF1,IF2,LDK,LDH,LDF
C 1 LPA,DM,DA,HELM,F1H,F1V,F2H,F2V,F,IFX,ISS,ISM,ISA,
C 2 T_A,T_RM,T_ALM,LDRA,LDR_RM,LDR_ALM,P_IH,C_IH,IER)
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-----------------------------------------------------------------------
SUBROUTINE T_CONV(EI,AKI,EF,AKF,QHKL,EN,HX,HY,HZ,IF1,IF2,
1 LDK,LDH,LDF,LPA,DM,DA,HELM,F1H,F1V,F2H,F2V,F,IFX,ISS,
2 ISM,ISA,T_A,T_RM,T_ALM,QM,LDRA,LDR_RM,LDR_ALM,P_IH,C_IH,IER)
C-----------------------------------------------------------------------
C
C INPUT
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 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 OUTPUT
C T_A TARGETS OF ANGLES A1-A6
C T_RM,T_ALM TARGETS OF RM ,LM
C QM TARGETS OF QM
C LDRA LOGICAL INDICATING WHICH ANGLES ARE TO BE DRIVEN
C LDR_RM,LDR_ALM LOGICAL INDICATING IF RM OR ALM ARE 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 DOUBLE PRECISION (A-H,O-Z)
PARAMETER(EPS1=1.D-1,EPS4=1.D-4)
C-----------------------------------------------------------------------
C PASSED PARAMETERS
REAL EI,AKI,EF,AKF,QHKL(3),EN,HX,HY,HZ
LOGICAL LDK(8),LDH,LDF,LPA
REAL DM,DA,HELM,F1H,F1V,F2H,F2V,F
REAL T_A(6),T_RM,T_ALM,QM
LOGICAL LDRA(6),LDR_RM,LDR_ALM
REAL P_IH(8),C_IH(4)
C-----------------------------------------------------------------------
C LOCAL VARIABLES
C
DIMENSION EDEF(2),AKDEF(2),DQHKL(3),DQT(3)
double precision ddm,dda
LOGICAL LMOAN(2),LQHKLE
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
IF(EDEF(1).LT.EPS1) GOTO 999
IER=0
AKDEF(1)=SQRT(EDEF(1)/F)
ELSE
IER=1
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
IF(EDEF(2).LT.EPS4) GOTO 999
IER=0
AKDEF(2)=SQRT(EDEF(2)/F)
ELSE
IER=1
IF(AKDEF(2).LT.EPS4) GOTO 999
IER=0
EDEF(2)=F*AKDEF(2)**2
ENDIF
EN=(3-2*IFX)*(EDEF(1)-EDEF(2))
ELSE IF (LQHKLE) THEN
LMOAN(3-IFX)=.TRUE.
EDEF(2)=EDEF(1)+(2*IFX-3)*EN
IER=1
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,IER)
IF (IER.EQ.0) THEN
AKI=DAKI
EI=DEI
T_A(1)=A1
T_A(2)=A2
T_RM=RM
T_ALM=ALM
LDRA(1)=.TRUE.
LDRA(2)=.TRUE.
LDR_RM=.TRUE.
LDR_ALM=.TRUE.
ELSE
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,IER)
IF (IER.EQ.0) THEN
AKF=DAKF
EF=DEF
T_A(5)=A5
T_A(6)=A6
LDRA(5)=.TRUE.
LDRA(6)=.TRUE.
ELSE
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
GOTO 999
ENDIF
ENDIF
C-----------------------------------------------------------------------
C DEAL WITH FLIPPERS AND HELMOTZ COILS IF LPA
C
IF (LPA.AND.(LMOAN(1).OR.LMOAN(2))) 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,
1 A4,QM,HELM,IER)
endif
C-----------------------------------------------------------------------
999 CONTINUE
IF (IER.NE.0) CALL ERRESO(IMOD,IER)
RETURN
END
SUBROUTINE EX_CASE(DX,ISX,AKX,AX1,AX2,RX,ALX,IER)
C-----------------------------------------------------------------------
C CALCULATE ANGLES ON MONO/ANALYSER
C CALCULATE AX1 AX2
C CALCULATE RX LX MONO CURVATURE AND LM FOR IN8
C
C INPUT
C DX D-SPACINGS
C ISX SENS OF SCATTERING ON CRYSTAL
C AKX TARGET OF MOMENTUM
C OUTPUT
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
C-----------------------------------------------------------------------
c Values of parameters
COMMON /CURVE/INX,C1RX,C2RX,RMIN,RMAX,CL1R
c
c INX=1 IN8 , INX=0 others instruments
C C1RX C2RX constants values to calculate RM on all instruments
C RMIN, RMAX min max on RNM
C CL1R constant value to calculate LM for IN8
c
C-----------------------------------------------------------------------
PARAMETER (PI=3.14159265358979323846264338327950E0)
PARAMETER(RD=57.29577951308232087679815481410517E0)
PARAMETER (EPS1=1.D-1)
C-----------------------------------------------------------------------
C PASSED PARAMETERS
DOUBLE PRECISION DX,AKX,AX1,AX2,ALX,RX
C-----------------------------------------------------------------------
C LOCAL VAR
DOUBLE PRECISION ARG,DC1RX,DC2RX,DRMIN,DRMAX,DCL1R
C-----------------------------------------------------------------------
C INIT AND TEST
C
IER=0
DC1RX=C1RX
DC2RX=C2RX
DRMIN=RMIN
DRMAX=RMAX
DCL1R=CL1R
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-----------------------------------------------------------------------
C CALCULATION OF THE TWO ANGLES
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 C1RX+C2RX/SIN(A1/RD)
C C1RX AND C2RX ARE CONSTANTS DEPENDING ON MONO AND MACHINES
C
C C1RX=.47
C C2RX=.244
C RMIN=0.
C RMAX=20.
C IN1/IN3/IN12/IN14/IN20 CASE
if (inx.EQ.0) then
RX=DMIN1(DMAX1(DC1RX+DC2RX/SIN(ABS(AX1)/RD),DRMIN),DRMAX)
else
C IN8 CASE
ALX=(DCL1R/SIN(AX2/RD))*COS(AX2/RD)
RX=DC2RX*SQRT(SIN(AX2/RD))-DC1RX
ENDIF
C-----------------------------------------------------------------------
999 CONTINUE
RETURN
END
C=========================================================================
SUBROUTINE SAM_CASE(QT,QM,QS,AKI,AKF,AX3,AX4,ISS,IER)
C-----------------------------------------------------------------------
C DEAL WITH SAMPLE ANGLES CALCULATION FROM Q VERTOR IN C-N PLANE
C CALCULATE A3 AND A4
C INPUT
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 OUTPUT
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 DOUBLE PRECISION (A-H,O-Z)
PARAMETER(RD=57.29577951308232087679815481410517E0)
PARAMETER (EPS3=1.D-3,EPS6=1.D-6)
C-----------------------------------------------------------------------
C PASSED PARAMETERS
DIMENSION QT(3)
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)/
1 (-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(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-----------------------------------------------------------------------
PARAMETER (PI=3.14159265358979323846264338327950E0)
PARAMETER(RD=57.29577951308232087679815481410517E0)
PARAMETER (EPS4=1.D-4)
C-----------------------------------------------------------------------
C PASSED PARAMETERS
DIMENSION T_IH(8),C_IH(4)
REAL*8 A4
C-----------------------------------------------------------------------
C INIT AND TEST
C
IER=1
IF (ABS(QM).LT.EPS4) goto 999
IER=0
DO IC=1,4
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
C HRAD !RADIAL COMP. OF H
C HDIR !DIRECTION OF H (IN RADIANS)
C HDIR2 !ANGLE BETWEEN FIELD AND AXE OF COIL 1
C
QPAR=AKI-AKF*COS(A4/RD)
QPERP=AKF*SIN(A4/RD)
PHI=ATAN2(QPAR,QPERP)
HRAD=SQRT(HX**2+HY**2)
IF(HRAD.GT.EPS4) HDIR=ATAN2(HY,HX)
HDIR2=PHI+HDIR+HELM/RD+0.5*PI
C-----------------------------------------------------------------------
C !CALC CURRENTS
C !POSITION OF PSP FOR COIL I
C
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)
C-----------------------------------------------------------------------
999 CONTINUE
RETURN
END
SUBROUTINE FLIP_CASE(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 PASSED PARAMETERS
DIMENSION T_IH(8)
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