*     *****************
      SUBROUTINE CALCUL
*     *****************
*
*     This subroutine calculates the structure factor and its
*     derivatives for reflection NQ.
*
      LOGICAL VERT
      REAL AS(3),AX(3),BS(3),BX(3),COSAR(30,24),DERIV(51),H(24,3),
     *HNN(3),RSA(30,3,3),RSAI(3,3),RSB(30,3,3),RSBI(3,3),SA(30),SB(30),
     *SINAR(30,24),SNEX(30),SNEXM(30),TEMP(30),TEXDER(15),TR(24),XI(14)
      INCLUDE 'COMMONS.FOR'
      COMMON/DERMAT/DERSTO(999,51)
      COMMON/DATA/CUT,NREFL,NH(3),IALFA,ISTAP,NCHAN,ICODE(999),
     *       HKL(3,999),SJJ(999),YOBS(15000),WEIGHT(15000),IORD(15000)
      COMMON/PARAMS/AINSTR(10),ACRYST(22), ATOM(14,30),
     *              LINSTR(10),LCRYST(22),LATOM(14,30),
     *              PINSTR(10),PCRYST(22),PATOM(14,30)
      COMMON/SPACE/JUMP,A,B,C,COSALF,COSBET,COSGAM,SINALF,SINBET,SINGAM
      SAVE
*
*     When angles in 1/10000 cycles, following statement should read
*     throughout this program:
*
*     DATA CONST/3183.09886184/ i.e. CONST = 10000/3.141592658...
      DATA CONST/11459.155903/
*                                Now CONST = 36000/3.141592658...
      DATA FACCUB/.66666666667/
*
      DATA GAUSS,TWOPI,SQRPI/0.9394372787,6.2831853072,1.7724585091/
*
*                        2 * sqrt(log 2/pi)   2 * pi       sqrt pi
*
      DATA SQRLN2/0.83255461116/,TR(1)/0.0/
*
      DO 500 NQ=1,NREFL
      DO 10 I=1,3
      HNN(I)=HKL(I,NQ)
   10 H(1,I)=HNN(I)
*
*     Calculation of q**2.
*
      QQ=PCRYST(1)*HNN(1)**2+PCRYST(2)*HNN(2)**2+PCRYST(3)*HNN(3)**2+
     *   PCRYST(4)*HNN(2)*HNN(3)+PCRYST(5)*HNN(1)*HNN(3)+
     *   PCRYST(6)*HNN(1)*HNN(2)
*
*     Overall temperature factor
*
      TAV=EXP(-PCRYST(7)*QQ/2)
*
*     Bragg's law: sin(theta)**2 = lambda**2 * q**2/4.
*
      SINTH=SLABDA*QQ
      COSTH=1-SINTH
      TANTH=SQRT(SINTH/COSTH)
*
*     BRAGG contains values of 2-theta in units specified by CONST.
*     GHALF = Gaussian halfwidth parameter.
*
      BRAGG(NQ)=CONST*ATAN(TANTH)
      GHALF=SQRT(PINSTR(3)/TANTH**2+TANTH*(PINSTR(4)*TANTH+PINSTR(5))+
     *           PINSTR(6))
*
*     Texture correction for each reflection.
*
      TEXCOR(NQ)=TEXTUR(NQ,TEXDER)
*
*     VERT = .TRUE. if asymmetry correction is to be calculated.
*
      VERT=BRAGG(NQ).LE.LIM
      DO 50 IR=2,IRL
      X=0
      DO 40 I=1,3
      Y=0
      X=EQUIV(4,I,IR)*HNN(I)+X
      DO 30 J=1,3
   30 Y=EQUIV(I,J,IR)*HNN(J)+Y
   40 H(IR,I)=Y
      TR(IR)=X
   50 CONTINUE
*
*     Calculation of COS(hx),SIN(hx) and temperature factor for
*     each atom.
*
      DO 90 I=1,N
      DO 60 J=1,14
   60 XI(J)=PATOM(J,I)
      DO 80 IR=1,IRL
      ARG=TR(IR)
      DO 70 J=1,3
   70 ARG=H(IR,J)*XI(J)+ARG
      ARG=TWOPI*ARG
      ARG2=H(IR,1)*H(IR,1)*XI( 9)+  H(IR,2)*H(IR,2)*XI(10)+
     1     H(IR,3)*H(IR,3)*XI(11)+2*H(IR,1)*H(IR,2)*XI(12)+
     2   2*H(IR,1)*H(IR,3)*XI(13)+2*H(IR,2)*H(IR,3)*XI(14)
      EXPARG=EXP(-ARG2)
      COSAR(IR,I)=COS(ARG)*EXPARG
   80 SINAR(IR,I)=SIN(ARG)*EXPARG
      TEMP(I)=EXP(-PATOM(4,I)*QQ/4)
   90 CONTINUE
*
*     Zeroize the derivatives of this reflection.
*
      DO 100 I=1,MAXS
  100 DERIV(I)=0
      IF(ICODE(NQ).EQ.2) GO TO 190
*
*     Calculate nuclear contribution to structure factor.
*
      AV=0
      BV=0
      AVIM=0
      BVIM=0
      DO 120 I=1,N
      IF(NTYP(I).EQ.0) GO TO 120
      SAI=0
      SBI=0
      DO 110 IR=1,IRL
      SAI=SAI+COSAR(IR,I)
      IF(ICENT.EQ.1) SBI=SINAR(IR,I)+SBI
  110 CONTINUE
      SAI=ICENT*SAI
      SA(I)=SAI
      SB(I)=SBI
      NI=NTYP(I)
*                                              X-rays
      IF(XRAYS.OR.ICODE(NQ).GE.4) THEN
         SNEXRE=(FF(NQ,NI)+DISPRE(NI))*PATOM(5,I)*TEMP(I)
         SNEXIM=DISPIM(NI)*PATOM(5,I)*TEMP(I)
         AVIM=SNEXIM*SAI+AVIM
         BVIM=SNEXIM*SBI+BVIM
         SNEXM(I)=2*SNEXIM*TAV
*                                              Neutrons
      ELSE
         SNEXRE=SCAT(NI)*PATOM(5,I)*TEMP(I)
         SNEXIM=0
      END IF
*
*     Calculate A and B of F
*
      AV=SNEXRE*SAI+AV
      BV=SNEXRE*SBI+BV
      SNEX(I)=2*SNEXRE*TAV
  120 CONTINUE
      FNN=(AV*AV+BV*BV+AVIM*AVIM+BVIM*BVIM)*TAV
      FNUC(NQ)=TEXCOR(NQ)*FNN
      IF(MAXS.EQ.0) GO TO 180
*
*     Nuclear contribution to derivatives - for atoms
*
      DO 160 I=1,N
      IF(NTYP(I).EQ.0) GO TO 160
      SNEXRE=SNEX(I)
      SNEXIM=SNEXM(I)
      SAI=SA(I)
      SBI=SB(I)
      DO 150 J=1,14
      K=LATOM(J,I)
      IF(K.EQ.0.OR.(J.GT.5.AND.J.LT.9)) GO TO 150
      IF(J.LE.3) THEN
*
*     Derivatives for atom coordinates - J = 1 to 3.
*
         SUMA=0
         SUMB=0
         DO 130 IR=1,IRL
         SUMA=H(IR,J)*SINAR(IR,I)+SUMA
         IF(ICENT.EQ.1) SUMB=H(IR,J)*COSAR(IR,I)+SUMB
  130    CONTINUE
         DER=-(SNEXRE*(AV*SUMA-BV*SUMB)+SNEXIM*(AVIM*SUMA-BVIM*SUMB))*
     *        ICENT*TWOPI
      ELSE IF(J.EQ.4) THEN
*
*     Derivatives for temperature factors  ---  J = 4.
*
         DER=-(SNEXRE*(AV*SAI+BV*SBI)+SNEXIM*(AVIM*SAI+BVIM*SBI))*QQ/4
      ELSE IF(J.EQ.5) THEN
*
*     Derivative for occupation number     ---  J = 5.
*
         DER=(SNEXRE*(AV*SAI+BV*SBI)+SNEXIM*(AVIM*SAI+BVIM*SBI))/
     *        PATOM(5,I)
      ELSE
*
*     Derivatives for anisotropic temperature parameters - J = 9 to 11.
*
         IF(J.LE.11) THEN
            K1=J-8
            K2=K1
         ELSE
            K1=J/14+1
            K2=J/13+2
         END IF
*
*     For the rest of the anisotropic temperature parameter -
*                                     J = 12 to 14.
*
         SUMA=0
         SUMB=0
         DO 140 IR=1,IRL
         SUMA=SUMA+H(IR,K1)*H(IR,K2)*COSAR(IR,I)
         IF(ICENT.EQ.1) SUMB=SUMB+H(IR,K1)*H(IR,K2)*SINAR(IR,I)
  140    CONTINUE
         DER=-(SNEXRE*(AV*SUMA+BV*SUMB)+SNEXIM*(AVIM*SUMA+BVIM*SUMB))*
     *        ICENT
         IF(K1.NE.K2) DER=2*DER
      END IF
      DERIV(K)=ATOM(J,I)*DER+DERIV(K)
  150 CONTINUE
  160 CONTINUE
*
*     Derivatives for - scale factor
*                     - overall temperature factor
*                     - texture parameters           - J = 8 to 22.
      K=LINSTR(2)
      IF(K.NE.0) DERIV(K)=DERIV(K)+AINSTR(2)*FNN*SCALE
      K=LCRYST(7)
      IF(K.NE.0) DERIV(K)=DERIV(K)-ACRYST(7)*FNN*QQ/2
      DO 170 J=8,22
      K=LCRYST(J)
      IF(K.NE.0) DERIV(K)=DERIV(K)+ACRYST(J)*FNN*TEXDER(J-7)
  170 CONTINUE
  180 IF(ICODE(NQ).EQ.1.OR.ICODE(NQ).GE.4) GO TO 450
*
*     Calculate magnetic contribution to structure factor and
*     derivatives.
*     For uniaxial spin configuration.
*
  190 IF(UNIAX) THEN
         X=FAC(NQ)
         FACCOS=(1+X)/2
         FACSIN=1-X
      END IF
      DO 200 I=1,3
      AX(I)=0
  200 BX(I)=0
      DO 260 I=1,N
      IF(MTYP(I).EQ.0) GO TO 260
      MEQI=MEQ(I)
      MTYPI=MTYP(I)
      SNEXRE=FF(NQ,MTYPI)*TEMP(I)*0.2695*PATOM(5,I)
      SNEX(I)=SNEXRE
      DO 220 J=1,3
      DO 220 K=1,3
      JK=(J-1)*3+K
      Y=0
      Z=0
      DO 210 IR=1,IRL
      X=RM(JK,MEQI,IR)
      Y=COSAR(IR,I)*X+Y
  210 IF(ICENT.EQ.1) Z=SINAR(IR,I)*X+Z
      RSAI(J,K)=Y
  220 RSBI(J,K)=Z
      DO 230 J=1,3
      DO 230 K=1,3
      RSAI(J,K)=RSAI(J,K)*ICENT
      RSA(I,J,K)=RSAI(J,K)
  230 RSB(I,J,K)=RSBI(J,K)
      DO 240 J=1,3
      AS(J)=0
      BS(J)=0
      DO 240 K=1,3
      AS(J)=PATOM(K+5,I)*RSAI(J,K)+AS(J)
      IF(ICENT.EQ.1) BS(J)=PATOM(K+5,I)*RSBI(J,K)+BS(J)
  240 CONTINUE
      DO 250 J=1,3
      AX(J)=AS(J)*SNEXRE+AX(J)
      IF(ICENT.EQ.1) BX(J)=BS(J)*SNEXRE+BX(J)
  250 CONTINUE
      SNEX(I)=SNEXRE*TAV
  260 CONTINUE
*
*     According to Halpern and Johnson
*
      IF(HALPER) THEN
         X=AX(1)**2+AX(2)**2+AX(3)**2+BX(1)**2+BX(2)**2+BX(3)**2+
     *     2*(COSALF*AX(2)*AX(3)+COSBET*AX(1)*AX(3)+COSGAM*AX(1)*AX(2)+
     *        COSALF*BX(2)*BX(3)+COSBET*BX(1)*BX(3)+COSGAM*BX(1)*BX(2))
         Y=(HNN(1)*AX(1)/A)**2+(HNN(2)*AX(2)/B)**2+(HNN(3)*AX(3)/C)**2+
     *     (HNN(1)*BX(1)/A)**2+(HNN(2)*BX(2)/B)**2+(HNN(3)*BX(3)/C)**2+
     *     2*HNN(1)*HNN(2)*(AX(1)*AX(2)+BX(1)*BX(2))/(A*B)+
     *     2*HNN(1)*HNN(3)*(AX(1)*AX(3)+BX(1)*BX(3))/(A*C)+
     *     2*HNN(2)*HNN(3)*(AX(2)*AX(3)+BX(2)*BX(3))/(B*C)
         FNN=X-Y/QQ
      END IF
*
*     For uniaxial spin configuration
*
      IF(UNIAX) FNN=FACSIN*(AX(3)**2+BX(3)**2)+
     *              FACCOS*(AX(1)**2+BX(1)**2+AX(2)**2+BX(2)**2+
     *                      2*COSGAM*(AX(1)*AX(2)+BX(1)*BX(2)))
*
*    For cubic spin configuration.
*
      IF(CUBIC) FNN=FACCUB*(AX(1)**2+AX(2)**2+AX(3)**2+
     *                      BX(1)**2+BX(2)**2+BX(3)**2)
*
*     Calculate magnetic structure factor.
*
      FNN=FNN*TAV
*
*     Total structure factor.
*
      FMAG(NQ)=TEXCOR(NQ)*FNN
      IF(MAXS.EQ.0) GO TO 500
*
*     Calculate derivatives.
*
      DO 300 I=1,3
      AS(I)=0
  300 BS(I)=0
      IF(HALPER) THEN
         X=2*(HNN(1)*AX(1)/A+HNN(2)*AX(2)/B+HNN(3)*AX(3)/C)/QQ
         AS(1)=2*(AX(1)+AX(2)*COSGAM+AX(3)*COSBET)-HNN(1)*X/A
         AS(2)=2*(AX(2)+AX(3)*COSALF+AX(1)*COSGAM)-HNN(2)*X/B
         AS(3)=2*(AX(3)+AX(1)*COSBET+AX(2)*COSALF)-HNN(3)*X/C
         IF(ICENT.EQ.1) THEN
            X=2*(HNN(1)*BX(1)/A+HNN(2)*BX(2)/B+HNN(3)*BX(3)/C)/QQ
            BS(1)=2*(BX(1)+BX(2)*COSGAM+BX(3)*COSBET)-HNN(1)*X/A
            BS(2)=2*(BX(2)+BX(3)*COSALF+BX(1)*COSGAM)-HNN(2)*X/B
            BS(3)=2*(BX(3)+BX(1)*COSBET+BX(2)*COSALF)-HNN(3)*X/C
         END IF
      END IF
      IF(UNIAX) THEN
         AS(1)=2*FACCOS*(AX(2)*COSGAM+AX(1))
         AS(2)=2*FACCOS*(AX(1)*COSGAM+AX(2))
         AS(3)=2*FACSIN*AX(3)
         IF(ICENT.EQ.1) THEN
            BS(1)=2*FACCOS*(BX(2)*COSGAM+BX(1))
            BS(2)=2*FACCOS*(BX(1)*COSGAM+BX(2))
            BS(3)=2*FACSIN*BX(3)
         END IF
      END IF
      IF(CUBIC) THEN
         DO 330 I=1,3
         AS(I)=2*AX(I)*FACCUB
         IF(ICENT.EQ.1) BS(I)=2*BX(I)*FACCUB
  330    CONTINUE
      END IF
*
*     Magnetic contribution to derivatives - for atoms.
*
      DO 430 I=1,N
      IF(MTYP(I).EQ.0) GO TO 430
      SNEXRE=SNEX(I)
      MEQI=MEQ(I)
      DO 420 J=1,8
      LPIJ=LATOM(J,I)
      IF(LPIJ.EQ.0) GO TO 420
      DO 340 K=1,3
      AX(K)=0
  340 BX(K)=0
*
*     Derivatives for atom coordinates - J = 1 to 3.
*
      IF(J.LE.3) THEN
         DO 370 L=1,3
         DO 360 K=1,3
         X=0
         Y=0
         DO 350 IR=1,IRL
         LM=(L-1)*3+K
         Z=RM(LM,MEQI,IR)
         X=H(IR,J)*SINAR(IR,I)*Z+X
         IF(ICENT.EQ.1) Y=H(IR,J)*COSAR(IR,I)*Z+Y
  350    CONTINUE
         AX(L)=PATOM(K+5,I)*X+AX(L)
         IF(ICENT.EQ.1) BX(L)=PATOM(K+5,I)*Y+BX(L)
  360    CONTINUE
         AX(L)=-AX(L)*SNEXRE*ICENT*TWOPI
         IF(ICENT.EQ.1) BX(L)=BX(L)*SNEXRE*TWOPI
  370    CONTINUE
*
*     Derivatives for temparature factor and occupation number -
*                                 J = 4 to 5.
      ELSE IF(J.LE.5) THEN
         DO 390 K=1,3
         DO 380 L=1,3
         AX(K)=RSA(I,K,L)*PATOM(L+5,I)+AX(K)
         IF(ICENT.EQ.1) BX(K)=RSB(I,K,L)*PATOM(L+5,I)+BX(K)
  380    CONTINUE
         AX(K)=AX(K)*SNEXRE*X1
         IF(ICENT.EQ.1) BX(K)=BX(K)*SNEXRE*X1
  390    CONTINUE
      ELSE
*
*     Derivatives for anisotropic temperature factors - J = 6 to 8.
*
         DO 400 K=1,3
         AX(K)=RSA(I,K,J-5)*SNEXRE
         IF(ICENT.EQ.1) BX(K)=RSB(I,K,J-5)*SNEXRE
  400    CONTINUE
      END IF
      DER=0
      DO 410 K=1,3
      IF(ICENT.EQ.1) THEN
         X1=BS(K)*BX(K)
      ELSE
         X1=0
      END IF
  410 DER=X1+AS(K)*AX(K)+DER
      DERIV(LPIJ)=ATOM(J,I)*DER+DERIV(LPIJ)
  420 CONTINUE
  430 CONTINUE
*
*     Derivatives for - scale factor
*                     - overall temperature factor
*                     - texture parameters           - J = 8 to 22.
      K=LINSTR(2)
      IF(K.NE.0) DERIV(K)=DERIV(K)+AINSTR(2)*FNN*SCALE
      K=LCRYST(7)
      IF(K.NE.0) DERIV(K)=DERIV(K)-ACRYST(7)*FNN*QQ/2
      DO 440 J=8,22
      K=LCRYST(J)
      IF(K.NE.0) DERIV(K)=DERIV(K)+ACRYST(J)*FNN*TEXDER(J-7)
  440 CONTINUE
  450 IF(MAXS.EQ.0) GO TO 500
*
*     Calculate derivatives for the rest of the parameters.
*     If FNN is less than 0.1 the derivatives can be approximated
*     to zero.
*
      FNN=FNUC(NQ)+FMAG(NQ)
      IF(FNN.LT.0.1) GO TO 480
*
*     Gaussian halfwidth parameter derivatives - first T, then U, V, W.
*
      X=TANTH*TANTH
      K=LINSTR(3)
      IF(K.NE.0) DERIV(K)=AINSTR(3)*FNN/(2*X*GHALF**2)+DERIV(K)
      DO 460 J=4,6
      K=LINSTR(J)
      IF(K.NE.0) DERIV(K)=X*AINSTR(J)*FNN/(2*GHALF**2)+DERIV(K)
  460 X=X/TANTH
*
*     Lorentzian halfwidth parameter derivatives.
*
      K=LINSTR(7)
      IF(K.NE.0) DERIV(K)=SQRLN2*TANTH*AINSTR(7)*FNN/GHALF+DERIV(K)
      K=LINSTR(8)
      IF(K.NE.0) DERIV(K)=SQRLN2*AINSTR(8)*FNN/(GHALF*SQRT(COSTH))+
     *                    DERIV(K)
*
*     Zeropoint derivatives and cell constant derivatives.
*
      K=LINSTR(9)
      IF(K.NE.0) DERIV(K)=AINSTR(9)*SQRPI*GAUSS*FNN/GHALF+DERIV(K)
      K=LINSTR(10)
      IF(K.NE.0) DERIV(K)=AINSTR(10)*SQRPI*GAUSS*FNN/(TANTH*GHALF)+
     *                    DERIV(K)
      DO 470 J=1,6
      K=LCRYST(J)
      IF(K.EQ.0) GO TO 470
      IF(J.LT.4) THEN
         X=HNN(J)*HNN(J)
      ELSE IF(J.EQ.4) THEN
         X=HNN(2)*HNN(3)
      ELSE IF(J.EQ.5) THEN
         X=HNN(1)*HNN(3)
      ELSE
         X=HNN(1)*HNN(2)
      END IF
      DERIV(K)=X*SLABDA*FNN*ACRYST(J)/(2*SQRT(SINTH*COSTH))+DERIV(K)
  470 CONTINUE
*
*     Asymmetry parameter derivative.
*
      K=LINSTR(1)
      IF((K.NE.0).AND.VERT) DERIV(K)=-AINSTR(1)*FNN/TANTH+DERIV(K)
*
*     Store derivatives for this reflection.
*
  480 DO 490 I=1,MAXS
  490 DERSTO(NQ,I)=DERIV(I)
  500 CONTINUE
      RETURN
      END
