             SUBROUTINE EXTICA(NONPOS,INEUTA,DIN,DOUT)
C
C     CALCULATION OF THE EXTINCTION IN  THE  MODEL OF  BECKER & COPPENS.
C     ALSO THE DERIVATIVES WORTH TO THE EXTINCTION PARAMETERS ARE  COM -
C                               PUTED.
C
C     WRITTEN BY UDO ZUCKER. ROUTINE TAKE PATTERN FROM XRAY 80 SYTEM AND
C                     FROM LINEX 74 PROGRAMME.
C
C
C     INPUT VARIABLES
C     ---------------
C
C     NONPOS = 1/2  FOR THE CURRENT REFLECTION PROCESSED THE EXTINCTION
C                   FACTOR   IS  NOT/ IS IMAGINARY .      (NON POSITIVE
C                   DEFINITE TENSOR OF EXTINCTION.)
C     INEUTA = 0/1  NEUTRON/XRAY RADIATION
C     DIN    =      DIRECTION COSINES OF DIRECT BEAM
C     DOUT   =      DIRECTION COSINES OF DIFFRACTED BEAM
C
C     IMPORTANT VARIABLES
C     -------------------
C
C     NIEXTA = TYPE OF EXTINCTION
C     IVERT = DISTRIBUTION FUNCTION SELECTOR
C     FCALC = IS HERE FCALC*FCALC
C
      DIMENSION VECTOR(3),TERM(3,3),DOUT(3),DIN(3),EXTVEC(3)
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'D.INC'
      INCLUDE 'E.INC'
      INCLUDE 'F.INC'
      INCLUDE 'G.INC'
      INCLUDE 'EXTN.INC'
      equivalence (gammap,exbeta)	
C
      EQUIVALENCE (QTQ(16),OMEGA1),(QTQ(17),OMEGA2)
     1,(QTQ(18),OMEGA3),(QTQ(12),ZP),(QTQ(13),ZS),(QTQ(14),ZPP)
     2,(QTQ(15),ZPS),(QTQ(10),DBRHO),(QTQ(11),DBG)
C
      NIEXTA = IBDF(1)+1
      IVERT = IBDF(2)
      IF(NIEXTA.GE.5)THEN
C
C---- CALCULATION OF UNIT VECTOR NORMAL TO THE PLANE DEFINED BY
C     VECTORS OF INCIDENT AND REFLECTED BEAM.
C
        EXTVEC(1)= DIN(2)*DOUT(3)-DIN(3)*DOUT(2)
        EXTVEC(2)=-DIN(1)*DOUT(3)+DIN(3)*DOUT(1)
        EXTVEC(3)= DIN(1)*DOUT(2)-DIN(2)*DOUT(1)
        DTEMP=1.0/SQRT(EXTVEC(1)**2+EXTVEC(2)**2+EXTVEC(3)**2)
        DO 300 J=1,3
          EXTVEC(J)=EXTVEC(J)*DTEMP
  300   CONTINUE
        DO 310 II=1,3
        IF(NIEXTA.EQ.6)THEN
          VECTOR(II)=DOUT(II)
        ELSE
          VECTOR(II)=EXTVEC(II)
        END IF
  310   CONTINUE
      END IF
C---- RHO = ( SIN(THETA)/LAMBDA )**2
      THETA=ASIN(SQRT(RHO)*WVLNTH)
      COS2T=COS(2*THETA)
      COS2T2=COS2T**2
      SIN2T=SQRT(1-COS2T2)
C
C---- GAMMS=GAMMA OF SECONDARY EXTINCTION  =  FUNCTION (TBAR)
C     (NOTATION LIKE THOSE OF THE BECKER AND COPPENS PUBLICATIONS.)
C---- EXY IS EXTINCTION FACTOR.
C
      GAMMS=EXBETA
      EXY=1.
C
C---- POLARISATION CORRECTION
C
      POL=HELPE*COS2T2
      POL=POL/(POL+HELPF)
      POLO=1.0-POL
C
C---- FUNCTIONS  A(THETA) ,  B(THETA) AND C OF PRIMARY EXTINCTION.
C
      AP=0.2+0.45*COS2T
      BP=0.22-0.12*(0.5-COS2T)**2
      CP=2.0
C---- SAME FOR SECONDARY EXTINCTION (G=GAUSSIAN, L=LORENTZIAN)
      AG=0.58+0.48*COS2T+0.24*COS2T2
      BG=0.02-0.025*COS2T
      AL=0.025+0.285*COS2T
      BL=0.15-0.2*(0.75-COS2T)**2
      IF(COS2T.LT.0.0)BL=-0.45*COS2T
C
C---- BRANCH CORRESPONDING TO EXTINCTION MODEL
C
      GO TO (660,400,390,400,320,320,320),NIEXTA
C
C---- ANISOTROPIC, CALCULATE VECTOR(TRANSPOSE)*TENSOR*VECTOR.
C
  320 DO 330 I=1,6
      IF(NIEXTA.EQ.6) THEN
        QTQ(I)=BDF(10+I)*WVLNTH**2
      ELSE
        QTQ(I)=BDF(4+I)
      END IF
  330 CONTINUE
      PROD=0.0
      DO 360 II=1,3
      DO 360 JJ=II,3
      IF(II.EQ.JJ)THEN
        CONST=1.0
        K=II
      ELSE
        CONST=2.0
        K=II+JJ+1
      END IF
      TERM(II,JJ)=VECTOR(II)*VECTOR(JJ)*CONST
      PROD=PROD+TERM(II,JJ)*QTQ(K)
  360 CONTINUE
C                 
      IF(PROD.GT.0.) GO TO 370
C---- RESULT IS NOT PHYSICAL.
      NONPOS=2
      HH=0.
      GO TO 380
C---- RESULT IS MEANINGFUL
  370 HH=1.0/SQRT(PROD)
  380 CONTINUE
      GO TO(580,400,390,400,400,390,400),NIEXTA
C---- COEFFICIENT FOR TYPE 2 EXTINCTION (ISOTROPIC - ANISOTROPIC)
  390 AS=AP
      BS=BP
      CS=2.0
      GO TO 410
C---- COEFFICIENTS FOR TYPE 1, DEPENDENT ON DISTRIBUTION TYPE.
  400 IF(IVERT.EQ.0) THEN
        AS=AG
        BS=BG
        CS=2.12
      ELSE
        AS=AL
        BS=BL
        CS=2.0
      END IF
C---- CALCULATE DERIVATIVES
  410 GO TO(580,430,420,440,430,420,440),NIEXTA
C---- TYPE 2  DERIVATIVE
  420 IF(NIEXTA.EQ.6) THEN
        RO=HH
      ELSE
        RO=BDF(4)/WVLNTH
      END IF
      BETA=RO
      DBRHO=1.0
      GO TO 500
C---- TYPE 1 DERIVATIVES DBRHO AND DBG
  430 IF(NIEXTA.EQ.5) THEN
        GPAR=HH
      ELSE
        GPAR=BDF(3)
      END IF
      BETA=GPAR/SIN2T
      DBG=1.0/SIN2T
      DBRHO=0.0
      RO=0.0
      GO TO 500
C---- GENERAL CASE
  440 RO=BDF(4)/WVLNTH
      IF(NIEXTA.EQ.4) THEN
        GPAR=BDF(3)
      ELSE
        GPAR=HH
      END IF
      IF(IVERT)460,460,480
C---- GENERAL CASE  DERIVATIVES (GAUSSIAN) DBRHO,DBG IN QTQ(10,11)
  460 IF(RO.LT.1.E-08.AND.GPAR.LT.1.E-08)GO TO 470
      BETA   = 1./SQRT(GPAR**2+(RO*SIN2T)**2)
      DBRHO= (BETA*GPAR)**3
      DBG= (BETA*RO)**3*SIN2T**2
      BETA   = GPAR*BETA*RO
      GO TO 500
C---- INITIALVALUES OF RO, GPAR (G,R) ARE ZERO.
  470 BETA=0.
      DBRHO=1.
      DBG=SIN2T**2
      GO TO 500
C---- SAME FOR LORENTZIAN DISTRIBUTION.
  480 IF(RO.LT.1.E-08.AND.GPAR.LT.1.E-08)GO TO 490
      BETA=1./(GPAR+RO*SIN2T)
      DBRHO=(BETA*GPAR)**2
      DBG=(BETA*RO)**2*SIN2T
      BETA=GPAR*BETA*RO
      GO TO 500
C---- INITIALVALUES ARE ZERO.
  490 BETA=0.
      DBRHO=1.
      DBG=SIN2T
  500 CONTINUE
C
C---- COMPUTE  XI (I=P,S) ; THERE IS NO POLARIZATION YET. START WITH
C     SECONDARY EXTINCTION.
      X2=GAMMS*BETA*FCALC
C---- PRIMARY  EXTINCTION, GAMMAP = CONSTANT COMPUTED IN SUBROUTINE RFINE.
      IF(NIEXTA.EQ.7.OR.NIEXTA.EQ.4) THEN
        XP=GAMMAP*FCALC*RO**2
        YP=ANA(AP,BP,CP,XP)
        ZP=ABLEIT(YP,AP,BP,CP,XP)
      ELSE
        XP=0.0
        YP=1.0
        ZP=0.0
      END IF
C---- NOW SECONDARY EXTINCTION. DERIVATIVE IN ZS
      XS=X2*YP
      YS=ANA(AS,BS,CS,XS)
      ZS=ABLEIT(YS,AS,BS,CS,XS)
C---- NOW POLARISATION EFFECTS IN CASE OF XRAY RADIATION.
      IF(INEUTA)520,510,520
C---- NEUTRON
  510 EXY=YP*YS
      OMEGA1=EXY*(1.+XP*ZP+XS*ZS*YP*(1.+XP*ZP))
      OMEGA2=EXY*YP*ZS
      OMEGA3=EXY*(2.*GAMMAP*RO*ZP+ZS*YP*(GAMMS*DBRHO+
     & 2.*X2*GAMMAP*RO*ZP))
      GO TO 530
C---- XRAY
  520 XPP=XP*COS2T2
C---- PRIMARY EXTINCTION
      YPP=ANA(AP,BP,CP,XPP)
C---- DERIVATIVES IN ZPP  (=0 FOR  TYPE 1 OR 2 )
      IF(NIEXTA.EQ.4.OR.NIEXTA.EQ.7) THEN
        ZPP=ABLEIT(YPP,AP,BP,CP,XPP)
      ELSE
        ZPP=0.0
      END IF
C---- SECONDARY EXTINC., DERIVATIVE IN ZPS
      XPS=YPP*X2*COS2T2
      YPS=ANA(AS,BS,CS,XPS)
      ZPS=ABLEIT(YPS,AS,BS,CS,XPS)
C---- COMPLETE EXTINCTION AND DERIVATIVES
      EXY=POLO*YP*YS+POL*YPP*YPS
      OMEGA1=POLO*YP*YS*(1.0+XP*ZP+XS*ZS*(1.0+XP*ZP))+
     &POL*COS2T2*YPP*YPS*((1.0/COS2T2)+XP*ZPP+YPP*X2*ZPS*(1.0+
     &XPP*ZPP))
      OMEGA2=POLO*YP*YS*ZS*YP+POL*COS2T2*YPP*YPS*ZPS*YPP
      OMEGA3=POLO*YP*YS*(2.0*GAMMAP*RO*ZP+ZS*YP*(GAMMS*
     &DBRHO+2.0*X2*GAMMAP*RO*ZP))+POL*COS2T2*YPP*YPS*(2.0*
     &GAMMAP*RO*ZPP+ZPS*YPP*(GAMMS*DBRHO+2.0*COS2T2*X2*
     &GAMMAP*RO*ZPP))
  530 CONTINUE
C---- GO AHEAD WITH XRAY AND NEUTRON. COMPUTE COMPLETE DERIVATIVES
C     MAXIMUM ARE SEVEN.
      IF(MREF.EQ.0) THEN
C---- REFINEMENT ON F
        QTQ(19)=0.5*((SQRT(FCALC))**3)/(SQRT(EXY))
      ELSE
C---- REFINEMENT ON F**2
        QTQ(19)=FCALC**2
      END IF
C
C---- NOW COMPUTATION IS SAME FOR BOTH CASES.
C
      GO TO(580,540,540,540,550,550,550),NIEXTA
C---- ISOTROPIC CASE ; DERIVATIVE IN QTQ(20) (+IN QTQ(21) IN G.CASE)
  540 IF(NIEXTA.NE.2) THEN
        QTQ(20)=QTQ(19)*OMEGA3/WVLNTH
      ELSE
        QTQ(20)=QTQ(19)*OMEGA2*GAMMS*DBG
      END IF
      IF(NIEXTA.EQ.4) QTQ(21)=QTQ(19)*OMEGA2*GAMMS*DBG
      GO TO 570
C---- ANISOTROPIC CASE. DERIVATIVES IN  QTQ(20 - 26)
  550 QTQ(27)=-0.5*HH/PROD
      DO 560 I=1,3
      DO 560 J=I,3
      IF(I.EQ.J) THEN
        K=I
      ELSE
        K=I+J+1
      END IF
      IF(NIEXTA.NE.6) THEN
        QTQ(19+K)=QTQ(19)*GAMMS*OMEGA2*DBG*TERM(I,J)*QTQ(27)
      ELSE
        QTQ(19+K)=QTQ(19)*OMEGA3*TERM(I,J)*QTQ(27)*WVLNTH**2
      END IF
C     QTQ(20)=QTQ(20)+QTQ(21)-0.5*QTQ(23)
  560 CONTINUE
      IF(NIEXTA.EQ.7) QTQ(26)=QTQ(19)*OMEGA3/WVLNTH
  570 CONTINUE
  580 CONTINUE
      IF(MREF.EQ.0) THEN
C---- REFINEMENT ON F
        VORFAK= OMEGA1/SQRT(FCALC*EXY)
      ELSE
C---- REFINEMENT ON F**2
        VORFAK=OMEGA1*2
      END IF
      IF(NIEXTA.LE.4) THEN
C---- ISOTROPIC CASES
        JUNK3 = 1
      ELSE
C---- ANISOTROPIC CASES
        JUNK3=6
      END IF
      IF(NIEXTA.EQ.4.OR.NIEXTA.EQ.7) JUNK3=JUNK3+1
      DO 650 II=1,JUNK3
  650 DRDE(II)=QTQ(19+II)
  660 RETURN
      END
      FUNCTION ANA(A,B,C,X)
C
C  COMPUTATION OF THE FUNCTION  YI(X)  (I=P,S) IN EXTINCTION
C   CORRECTION OF BECKER AND COPPENS
      ANA=1.0/SQRT(1.0+C*X+A*X**2/(1.0+B*X))
      RETURN
      END
C          
      FUNCTION ABLEIT(YI,A,B,C,X)
C
C COMPUTE   1/YI(X)*DERIVATIVE(YI(X)) WTH. TO X   (I=P,S)
C
      ABLEIT=-0.5*YI**2*(C+A*X*(2.0+B*X)/(1.0+B*X)**2)
      RETURN
      END
