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