             SUBROUTINE EXCALC(NONPOS,HELPE,HELPF)
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     HELPE  =      COEFFICIENT OF POLARISATION CORRECTION,  CALCULATED
C                   IN ROUTINE REFINE.
C     HELPF  =      SAME AS HELPE
C
C     IMPORTANT VARIABLES
C     -------------------
C
C     NIEXT = TYPE OF EXTINCTION
C     NDISTR = DISTRIBUTION FUNCTION SELECTOR
C     FCALC = IS HERE FCALC*FCALC
C
C
      DIMENSION VECTOR(3),TERM(3,3),QTQ(35),IQTQ(35)
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'D.INC'
      INCLUDE 'E.INC'
      INCLUDE 'G.INC'
      INCLUDE 'F.INC'
      INCLUDE 'EXTN.INC'
C
      EQUIVALENCE (QTQ,IQTQ)
C
      NONPOS=1
      IF(NIEXT.NE.0)GO TO 9000
C
C ZACHARIASEN ISOTROPIC EXTINCTION HERE
C
      EXX=EXBETA*EXTIN(1)*FCALC
      DRDE(1)=1.0
      RETURN
C
C BECKER-COPPENS FORMULATION
C
C---- THE VALUE GAMMAS IS IN BINARY STORAGE CALLED EXBETA.  EXX IS EX -
C     TINCTION FACTOR.
C
9000  GAMMS=EXBETA
      EXX=1.
C
C---- COEFFICIENT TO FIT THE VALUE TO THE PROGRAMME SYSTEM.
C
      VORFAK=1./SQRT(FCALC)
C---- RHO = ( SIN(THETA)/LAMBDA )**2
      THETA=ASIN(SQRT(RHO)*WVLN)
      COS2T=COS(2*THETA)
      COS2T2=COS2T**2
      SIN2T=SQRT(1-COS2T2)
C
C---- POLARISATION CORRECTION
C
      POL=HELPE*COS2T2
      POL=POL/(POL+HELPF)
      POLO=1.0-POL
C
C---- FUNCTIONS  A(THETA) ,  B(THETA) UND 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
      GOTO(5719,5716,5719,5690,5690,5690),NIEXT
C
C---- ANISOTROPIC, CALCULATE VECTOR(TRANSPOSE)*TENSOR*VECTOR.
C
 5690 DO 5695 I=1,6
 5695 QTQ(I)=EXTIN(I)
      IF(NIEXT.NE.5) GOTO 5703
      DO 5700 I=1,6
 5700 QTQ(I)=BDF(10+I)*WVLN**2
 5703 DO 5706 II=1,3
      VECTOR(II)=EXANG(3+II)
      IF(NIEXT.EQ.5)VECTOR(II)=EXANG(II)
 5706 CONTINUE
      PROD=0.0
      DO 5709 II=1,3
      DO 5709 JJ=II,3
      CONST=1.0
      IF(II.NE.JJ)CONST=2.0
      K=II+JJ+1
      IF(II.EQ.JJ)K=II
      TERM(II,JJ)=VECTOR(II)*VECTOR(JJ)*CONST
      PROD=PROD+TERM(II,JJ)*QTQ(K)
 5709 CONTINUE
C
      IF(PROD.LE.0.) THEN
C
C---- RESULT HAS NO PHYSICAL MEANING
        NONPOS=2
        H=0.
      ELSE
C
C---- RESULT IS MEANINGFUL
C
        H=1.0/SQRT(PROD)
      END IF
C
C---- HERE NIEXT IS GREATER 3
C
      GOTO(5719,5716,5719,5719,5716,5719),NIEXT
C
C---- COEFFICIENT FOR TYPE 2 EXTINCTION (ISOTROPIC - ANISOTROPIC)
C
 5716 AS=AP
      BS=BP
      CS=2.0
      GOTO 5723
C
C---- COEFFICIENTS FOR TYPE 1, DEPENDENT ON DISTRIBUTION TYPE.
C
 5719 IF(NDISTR.EQ.0)THEN
        AS=AG
        BS=BG
        CS=2.12
      ELSE
        AS=AL
        BS=BL
        CS=2.0
      END IF
C
C---- CALCULATE DERIVATIVES
C
5723  GO TO(5729,5726,5736,5729,5726,5736),NIEXT
C
C---- TYPE 2  DERIVATIVE
C
 5726 RO=BDF(4)/WVLN
      IF(NIEXT.EQ.5) RO=H
      BETA=RO
      QTQ(10)=1.0
      GOTO 5749
C
C---- TYPE 1 DERIVATIVES DBRHO AND DBG IN QTQ(10,11)
C
 5729 GPAR=BDF(3)
      IF(NIEXT.EQ.4) GPAR=H
      BETA=GPAR/SIN2T
      QTQ(11)=1.0/SIN2T
      QTQ(10)=0.0
      RO=0.0
      GOTO 5749
C
C---- GENERAL CASE
C
 5736 RO=BDF(4)/WVLN
      GPAR=BDF(3)
      IF(NIEXT.EQ.3) GOTO 5739
      GPAR=H
 5739 IF(IVERT)5743,5743,5746
C
C---- GENERAL CASE  DERIVATIVES (GAUSSIAN) DBRHO,DBG IN QTQ(10,11)
C
 5743 IF(RO.LT.1.E-08.AND.GPAR.LT.1.E-08)GOTO 5747
      BETA   = 1./SQRT(GPAR**2+(RO*SIN2T)**2)
      QTQ(10)= (BETA*GPAR)**3
      QTQ(11)= (BETA*RO)**3*SIN2T**2
      BETA   = GPAR*BETA*RO
      GOTO 5749
C
C---- INITIALVALUES OF RO, GPAR (G,R) ARE ZERO.
C
 5747 BETA=0.
      QTQ(10)=1.
      QTQ(11)=SIN2T**2
      GOTO 5749
C
C---- SAME FOR LORENTZIAN DISTRIBUTION.
C
 5746 IF(RO.LT.1.E-08.AND.GPAR.LT.1.E-08)GOTO 5748
      BETA=1./(GPAR+RO*SIN2T)
      QTQ(10)=(BETA*GPAR)**2
      QTQ(11)=(BETA*RO)**2*SIN2T
      BETA=GPAR*BETA*RO
      GOTO 5749
C---- INITIALVALUES ARE ZERO.
 5748 BETA=0.
      QTQ(10)=1.
      QTQ(11)=SIN2T
 5749 CONTINUE
C
C---- COMPUTE  XI (I=P,S) ; THERE ARE NO POLARISATION YET. START WITH
C     SECONDARY EXTINCTION.
C
      X2=GAMMS*BETA*FCALC
      X1=0.0
C
C---- PRIMARY  &  SEKUNDARY  EXT., GAMMAP = CONSTANT COMPUTED IN SUB-
C     ROUTINE REFINE.
C
      IF(NIEXT.EQ.6.OR.NIEXT.EQ.3)     X1=GAMMAP*FCALC*RO**2
      XP=X1
C
C---- AMOUNT OF PRIMARY EXTINCTION TO COMPLETE EXTINCTION. DERIVATIVE
C     IN QTQ(12)
C
      YP=ANA(AP,BP,CP,XP)
C
C---- QTQ(12)=0  IF NOT GENERAL CASE.
C
      IF(NIEXT.EQ.3.OR.NIEXT.EQ.6) THEN
        QTQ(12)=ABLEIT(AP,BP,CP,XP)
      ELSE
        QTQ(12)=0.0
      END IF
C
C---- NOW SECONDARY EXTINCTION. DERIVATIVE IN QTQ(13)
C
      XS=X2*YP
      YS=ANA(AS,BS,CS,XS)
      QTQ(13)=ABLEIT(AS,BS,CS,XS)
C
C---- NOW POLARISATION EFFECTS IN CASE OF XRAY RADIATION.
C
      IF(INEUT)5756,5753,5756
C
C---- NEUTRON
C
 5753 EXX=YP*YS
      QTQ(16)=EXX*(1.+XP*QTQ(12)+XS*QTQ(13)*(1.+XP*QTQ(12)))
      QTQ(17)=EXX*YP*QTQ(13)
      QTQ(18)=EXX*(2.*GAMMAP*RO*QTQ(12)+QTQ(13)*YP*(GAMMS*QTQ(10)+
     & 2.*X2*GAMMAP*RO*QTQ(12)))
      GOTO 5759
C
C---- XRAY
C
 5756 XPP=X1*COS2T2
C
C---- PRIMARY EXTINCTION
C
      YPP=ANA(AP,BP,CP,XPP)
C---- DERIVATIVES IN QTQ(14)  (=0 FOR  TYPE 1 OR 2 )
      IF(NIEXT.EQ.3.OR.NIEXT.EQ.6) THEN
        QTQ(14)=ABLEIT(AP,BP,CP,XPP)
      ELSE
        QTQ(14)=0.0
      END IF
C
C---- SECONDARY EXTINC., DERIVATIVE IN QTQ(15)
C
      XPS=YPP*X2*COS2T2
      YPS=ANA(AS,BS,CS,XPS)
      QTQ(15)=ABLEIT(AS,BS,CS,XPS)
C
C---- COMPLETE EXTINCTION AND DERIVATIVES
C
      EXX=POLO*YP*YS+POL*YPP*YPS
C
C---- QTQ(16)=QMEGA1 , QTQ(17)=OMEGA2 , QTQ(18)=OMEGA3   (NOTATION OF
C     BECKER AND COPPENS PUBLICATION)
C
      QTQ(16)=POLO*YP*YS*(1.0+XP*QTQ(12)+XS*QTQ(13)*(1.0+XP*QTQ(12)))+
     &POL*COS2T2*YPP*YPS*((1.0/COS2T2)+XP*QTQ(14)+YPP*X2*QTQ(15)*(1.0+
     &XPP*QTQ(14)))
      QTQ(17)=POLO*YP*YS*QTQ(13)*YP+POL*COS2T2*YPP*YPS*QTQ(15)*YPP
      QTQ(18)=POLO*YP*YS*(2.0*GAMMAP*RO*QTQ(12)+QTQ(13)*YP*(GAMMS*
     &QTQ(10)+2.0*X2*GAMMAP*RO*QTQ(12)))+POL*COS2T2*YPP*YPS*(2.0*
     &GAMMAP*RO*QTQ(14)+QTQ(15)*YPP*(GAMMS*QTQ(10)+2.0*COS2T2*X2*
     &GAMMAP*RO*QTQ(14)))
C
 5759 CONTINUE
C
C---- GO AHEAD WITH XRAY AND NEUTRON. COMPUTE COMPLETE DERIVATIVES
C     MAXIMUM ARE SEVEN.
C
      IF(MREF.EQ.0) THEN
C---- REFINEMENT ON F
        QTQ(19)=0.5*((SQRT(FCALC))**3)/(SQRT(EXX))
      ELSE
C---- REFINEMENT ON F**2
        QTQ(19)=FCALC**2
      END IF
C
C---- NOW COMPUTATION IS SAME FOR BOTH CASES.
C
      GO TO (5763,5763,5763,5766,5766,5766),NIEXT
C
C---- ISOTROPIC CASE ; DERIVATIVE IN QTQ(20) (+IN QTQ(21) IN G.CASE)
C
 5763 IF(NIEXT.NE.1)THEN
        QTQ(20)=QTQ(19)*QTQ(18)/WVLN
      ELSE
        QTQ(20)=QTQ(19)*QTQ(17)*GAMMS*QTQ(11)
      END IF
      IF(NIEXT.EQ.3) QTQ(21)=QTQ(19)*QTQ(17)*GAMMS*QTQ(11)
      GOTO 5773
C
C---- ANISOTROPIC CASE. DERIVATIVES IN  QTQ(20 - 26)
C
 5766 QTQ(27)=-0.5*H/PROD
      DO 5769 I=1,3
      DO 5769 J=I,3
      K=I+J+1
      IF(I.EQ.J) K=I
      IF(NIEXT.NE.5)THEN
        QTQ(19+K)=QTQ(19)*GAMMS*QTQ(17)*QTQ(11)*TERM(I,J)*QTQ(27)
      ELSE
        QTQ(19+K)=QTQ(19)*QTQ(18)*TERM(I,J)*QTQ(27)*WVLN**2
      END IF
      QTQ(20)=QTQ(20)+QTQ(21)-0.5*QTQ(23)
 5769 CONTINUE
      IF(NIEXT.EQ.6) QTQ(26)=QTQ(19)*QTQ(18)/WVLN
 5773 CONTINUE
 5901 CONTINUE
C
C---- NOW COMPUTE MODIFICATIONS OF THE DERIVATIVES WTH TO OTHER PARA -
C     METERS. DISTINGUISH BETWEEN REFINEMENT ON F AND ON F*F.
C
      IF(MREF.EQ.0) THEN
C---- REFINEMENT ON F   >>>> VORFAKTOR=  OMEGA1 / FCALC*SQRT(EXX)
        VORFAK= QTQ(16)/SQRT(FCALC*EXX)
      ELSE
C *** REFINEMENT ON F**2 >>>> VORFAKTOR= OMEGA1*2
        VORFAK=QTQ(16)*2
      END IF
C
C---- CLEAR SYSTEM ARRAY FOR DERIVATIVES WTH TO EXTINCTION PARAMETRS.
C
      DO 6000 II=1,7
 6000 DRDE(II)=0.0
C
C---  STORE DERIVATIVES IN ACCORDANCE TO PROMETHEUS DESCRIBTION.
C
      GOTO (5956,5956,5956,5963,5963,5963),NIEXT
C
C---- ISOTROPIC CASES
C
 5956 JUNK3 = 1
      IF(NIEXT.EQ.3) JUNK3=2
      DO 5959 II=1,JUNK3
 5959 DRDE(II)=QTQ(19+II)
      GOTO 5999
C
C---- ANISOTROPIC CASES
C
 5963 JUNK3=6
      IF(NIEXT.EQ.6) JUNK3=7
      DO 5966 II=1,JUNK3
 5966 DRDE(II)=QTQ(19+II)
 5999 RETURN
      END
      FUNCTION ANA(A,B,C,X)
C
C     COMPUTATION OF THE FUNCTION  YI(X)  (I=P,S) IN EXTINCTIONCOR-
C     RECTION OF BECKER U. COPPENS
      ANA=1.0/SQRT(1.0+C*X+A*X**2/(1.0+B*X))
      RETURN
      END
      FUNCTION ABLEIT(A,B,C,X)
C
C     COMPUTE   1/YI(X)*DERIVATIVE(YI(X)) WTH. TO X   (I=P,S)
C
      ABLEIT=-0.5*ANA(A,B,C,X)**2*(C+A*X*(2.0+B*X)/(1.0+B*X)**2)
      RETURN
      END
