410 lines
13 KiB
Fortran
Executable File
410 lines
13 KiB
Fortran
Executable File
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
|