         SUBROUTINE EXTINK(RNEUEX,NIEXT,NDISTR,KCYL)
C
C      READING OF PARAMETERS OF THE EXTINCTION MODEL  OF BECKER & COPPENS
C      AS WELL AS CONVERTING AND PRINTING OF VALUES. PROGRAMME TAKES PAT-
C                    TERN OF THE PROGRAMME XRAY 80.
C
C                         WRITTEN BY UDO ZUCKER
C
C      INPUT VARIABLES
C      ---------------
C
C      RNEUEX = CARRY INFORMATION ABOUT THE EXTINCTION PARAMETERS JUST READ
C      NIEXT  = NEW EXTINCTION TYPE.
C      NDISTR = NEW DISTRIBUTION FUNCTION.
C      KCYL   = NUMBER OF CURRENT CYCLE.
C
      DIMENSION RNEUEX(7)
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'EXTN.INC'
      INCLUDE 'IALPH.INC'
      INCLUDE 'E.INC'
C
   10 FORMAT('0EXTINCTION PARAMETERS')
   12 FORMAT(13H0ERRORMESSAGE/1X,12(1H-)//    22H ONE OR MORE NEGATIVE ,
     & 58HISOTROPIC EXTINCTION PARAMETERS WERE ENCOUNTERED ON INPUT./
     & /32H THESE HAVE BEEN RESET POSITIVE.//1H )
   16 FORMAT(8H0WARNING/1X,7(1H-)//32H OLD AND NEW EXTINCTIONTYPES DO ,
     & 15HNOT CORRESPOND./42H0NO CONVERSION POSSIBLE. NEW TYPE USED AS ,
     & 13HINITIAL TYPE.//1H )
   18 FORMAT(13H0ERRORMESSAGE/1X,12(1H-)//  22H ISOTROPIC EXTINCTION ,
     & 58HPARAMETER IS ZERO. NO CONVERSION POSSIBLE. NO EXTINCTION ,
     & 11HCALCULATED.)
   20 FORMAT(1H0,19A4)
   22 FORMAT(1H0,34X,10H    R     ,2X,10H Z11 OR G ,2X,10H    Z22   ,2X,
     &10H    Z33   ,2X,10H    Z12   ,2X,10H    Z13   ,2X,10H    Z23   )
   24 FORMAT(1H0,8A4,12X,F10.4)
   26 FORMAT(1H0,8A4,F10.4)
   28 FORMAT(1H0,8A4,F10.4,2X,F10.4)
   30 FORMAT(1H0,8A4,12X,6(2X,E10.4))
   32 FORMAT(1H0,8A4,F10.4,2X,6(1X,E10.4))
   34 FORMAT(1H0,8A4,49X,E10.4)
   36 FORMAT(1H0,8A4,12X,3(2X,E10.4))
   38 FORMAT(1H0,8A4,2X,E10.4)
   42 FORMAT('0EIGENVALUES OF EXTINCTION TENSOR',6X,3E14.4)
   44 FORMAT(1H0,8A4,6X,3E14.4)
C
C---- PRINT HEADING.
C
      WRITE(IOUT,10)
C---- INITIALIZE
      PI   = ACOS(-1.)
      JUNK5= 0
C
C---- REDUCE VALUES OF LEGAL RANGE.
C
      NIEXT=MAX0(MIN0(NIEXT,6),0)
      NDISTR=MAX0(MIN0(NDISTR,1),0)
C
C---- CLEAR ARRAY FOR COMPUTATION.
      DO 300 I=1,16
  300 QTQ(I)=0.0
C
C---- ARE THE PARAMETERS READ IN ZERO, CONVERSION OF EXTINCTION TYPES IS SUPPOSED.
C
      IF(RNEUEX(1).LT.1.E-06) GOTO 400
C
C---- FIRST CYCLE.
C
      IQTQ(1)=NIEXT
      IQTQ(2)=NDISTR
      NIEXTA  =NIEXT+1
C
  310 GO TO (630,320,320,330,360,360,360),NIEXTA
C
C---- ISOTROPIC EXTINCTION
C
  320 IF(NIEXTA.EQ.3) THEN
        QTQ(4)=ABS(RNEUEX(1))
      ELSE
        QTQ(3)=ABS(RNEUEX(1))
      END IF
      GO TO 340
  330 QTQ(3)=ABS(RNEUEX(1))
      QTQ(4)=ABS(RNEUEX(2))
C
C---- WARNING, IF INPUT PARAMETER WAS NEGATIVE.
C
  340 IF(RNEUEX(1).LT.0.0.OR.RNEUEX(2).LT.0.0) GOTO 350
      GOTO 630
C
C---- VALUE IS NEGATIVE.
C
  350 WRITE(IOUT,12)
      GOTO 630
C
C---- ANISOTROPIC EXTINCTION.
C
  360 DO 370 I=1,6
      IF(NIEXTA.NE.6) THEN
        QTQ(I+4)  = RNEUEX(I)
      ELSE
        QTQ(I+10) = RNEUEX(I)
      END IF
  370 CONTINUE
      IF(NIEXTA.EQ.7) QTQ(4)    = ABS(RNEUEX(7))
C
C---- CHECK IN GENERAL CASE R FOR POSITIVITY.
C
  390 IF(NIEXTA.NE.7.OR.RNEUEX(7).GT.0.0) GOTO 630
C
C---- NEGATIVE VALUE FOUND.
C
      WRITE(IOUT,12)
      GOTO 630
C
C---- TO THIS POINT BRANCH LED, IF NO NEW EXTINCTION  PARAMETERS OF
C     VALUES DIFFERENT FROM ZERO WERE READ. CONVERSION IS SELECTED.
C
  400 CONTINUE
C
C---- READ  DATAS OF LAST CYCLE.
C
      IQTQ(1) = IBDF(1)
      NIEXTA   = IQTQ(1)+1
      IQTQ(2) = IBDF(2)
C
      GOTO(440,410,410,410,420,420,420),NIEXTA
C
C---- ISOTROPIC EXTINCTION.
C
  410 IF(NIEXTA.NE.3)QTQ(3)=BDF(3)
      IF(NIEXTA.NE.2)QTQ(4)=BDF(4)
      GOTO 440
C
C---- ANISOTROPIC EXTINCTION.
C
  420 DO 430 I=1,6
      IF(NIEXTA.NE.6) QTQ(I+4)  = BDF(I+4)
      IF(NIEXTA.EQ.6) QTQ(I+10) = BDF(I+10)
  430 CONTINUE
      IF(NIEXTA.EQ.7) QTQ(4) = BDF(4)
  440 CONTINUE
C
C---- CONVERSION.
C
  450 IF(IQTQ(1).NE.NIEXT)GOTO 460
C
C---- ONLY DISTRIBUTION FUNCTION TO CHANGE.
C
      IQTQ(2)=NDISTR
      GOTO 630
C
C---- NO EXTINCTION CORRECTION SELECTED , PARAMETER = 0
C
  460 IF(NIEXT.EQ.0) GOTO 610
C
C---- CHECK FOR CORRESPONDENCE.
C
      IF(IQTQ(1).EQ.0)GOTO 470
      IF(IABS(IQTQ(1)-NIEXT).EQ.3) GOTO 480
C
C---- THERE IS NO CONVERSION POSSIBLE, BECAUSE OF NONCORRESPONDENCE.
  470 WRITE(IOUT,16)
C
C---- RESET ANISOTROPIC TO ISOTROPIC AND TAKE ISOTROPIC MODEL
C
      IF(NIEXT.GT.3)NIEXT=NIEXT-3
      NIEXTA=NIEXT+1
      IQTQ(1)=NIEXT
      IQTQ(2)=NDISTR
      GOTO 310
C
  480 ASSIGN 590 TO NTRET
      ASSIGN 630 TO NRETN
C
C---- ASSIGN ADRESSES FOR BRANCHES.
C
      IF(NIEXT.LE.3) ASSIGN 520 TO NTRET
      IF(NIEXT.GE.4) ASSIGN 490 TO NTRET
      JUNK3=3
      JUNK4=4
      IF(NIEXT.EQ.2) THEN
        JUNK3=4
        JUNK4=10
      END IF
      IF(NIEXT.EQ.5) THEN
        JUNK3=4
        JUNK4=10
      END IF
C
C---- GO TO NTRET.
C
      GOTO NTRET,(490,520,590)
C
C---- CONVERSION   ISOTROPIC TO ANISOTROPIC
C
  490 IF(QTQ(JUNK3).GT.0.0) GOTO 500
C
C---- BUT ONLY, IF PARAMETER IS NOT ZERO.
C
      ASSIGN 600 TO NRETN
      GOTO 590
C
C---- YES IT IS, CALCULATE EQUIVALENT TENSOR.
C
  500 QTQ(JUNK3)= 1.0/QTQ(JUNK3)
      DO 510 I=1,3
  510 QTQ(JUNK4+I)=QTQ(JUNK3)**2
      QTQ(JUNK4+4)=0.0
      QTQ(JUNK4+5)=0.0
      QTQ(JUNK4+6)=0.0
C---- DELETE OLD INFORMATION
C
      QTQ(JUNK3)=0.0
      ASSIGN 630 TO NRETN
      GOTO 590
C
C---- ANISOTROPIC EXT. >>>> ISOTROPIC
C
  520 CONTINUE
C
C---- TRANSFER  QTQ(3 - 16) IN QTQ(40 - 53)
C
      DO 530 I=1,14
  530 QTQ(I+39)=QTQ(I+2)
C
C---- THIS HAD BEEN DONE BECAUSE JACOBI SHOULD BE READ.
C
      JUNK4=JUNK4+37
C
C---- ORTHOGONAL EXTINCTION TENSOR. RESULT IN QTQ(1)-QTQ(9)
C
  540 DO 550 I=1,3
      DO 550 J=1,3
      K=3*(I-1)+J
      IF(I.EQ.J)THEN
        QTQ(K)=QTQ(JUNK4+I)
      ELSE
        QTQ(K)=QTQ(JUNK4+1+I+J)
      END IF
  550 CONTINUE
C
C *** CALL JACOBI-SUBROUTINE TO CALCULATE EIGENVALUES (QTQ(1,5,9))
C
      CALL JACOBI(QTQ)
  560 IF(JUNK5.EQ.1) GOTO 890
      QTQ(54)=SQRT((QTQ(1)+QTQ(5)+QTQ(9))/3.0)
C
C---- ISOTROPIC EXTINCT. VALUE IN QTQ(54)
C
      IF(QTQ(54).LT.1.E-08) GOTO 570
      QTQ(54)=1.0/QTQ(54)
  570 JUNK4=JUNK4-37
C
C---- NEW PARAM. AGAIN  IN QTQ(JUNK3)
C
      QTQ(JUNK3)=QTQ(54)
      IF(NIEXT.EQ.3) QTQ(4)=QTQ(41)
      DO 580 I=1,6
  580 QTQ(JUNK4+I)=0.0
C
C---- FINISH OF CONVERSION
      ASSIGN 630 TO NRETN
  590 CONTINUE
C
C---- FILL IN IQTQ(1 - 2)  WITH NEW VALUES.
C
      IQTQ(1)=NIEXT
      IQTQ(2)=NDISTR
C
C---- IF ALL CORRECT GO ON WITH PRINTING.
C
      GOTO NRETN,(600,630)
C
C---- THERE WAS AN ERROR.
C
  600 CONTINUE
      WRITE(IOUT,18)
C
C---- RESET TO NO EXTINCTION
C
  610 DO 620 I=3,16
  620 QTQ(I)=0.0
      IQTQ(1)=0
      IQTQ(2)=0
  630 CONTINUE
C
C---- NOW FILL IN NEW VALUES IN FIELD BDF, STORED IN MAIN FILE.
C
      IBDF(1) = IQTQ(1)
      IBDF(2) = IQTQ(2)
      DO 640 I=3,16
      BDF(I)   =QTQ(I)
  640 QTQ(I+37)=QTQ(I)
C
C---- EASIER NAME FOR PRINTING SECTION.
C
      ITYP  =IQTQ(1)+1
C
C
C---- PRINT SECTION ------------------------------------------------
C
C---- CLEAR FIELD WITH BLANKS
C
      DO 650 I=1,19
  650 IQTQ(I)=IDATA1(27)
      IF(ITYP.NE.1)GOTO 670
C---- FILL IN *NO EXTINCTION CORRECTION*
      DO 660 I=1,7
  660 IQTQ(I)=IDATA1(I)
      GOTO 800
  670 CONTINUE
C---- FILL IN *ANISOTROPIC*
      DO 680 I=1,3
  680 IQTQ(I)=IDATA1(I+8)
C---- CHANGE TO *ISOTROPIC*, IF ISOTROPIC CASE COMPUTED.
      IF(ITYP.LE.4)IQTQ(1)=IDATA1(8)
C---- FILL IN *EXTINCTION CORRECTION*
      DO 690 I=1,6
  690 IQTQ(I+3)=IDATA1(I+1)
C
      GOTO(800,700,740,750,700,740,750),ITYP
C
C---- FILL IN *TYPE 1 GAUSSIAN DISTRIBUTION*
  700 DO 710 I=1,2
      IQTQ(I+9)=IDATA1(I+11)
      IQTQ(I+11)=IDATA1(I+14)
      IQTQ(I+13)=IDATA1(I+19)
  710 IQTQ(I+15)=IDATA1(I+21)
C---- IF LORENTZIAN COMPUTED, CHANGE GAUSSIAN TO LORENTZIAN
      IF(NDISTR)800,800,720
  720 DO 730 I=1,7
  730 IQTQ(I+11)=IDATA1(I+16)
      GOTO 800
C---- FILL IN *ISOTROPIC TYPE 2 EXTINCTION*
  740 IQTQ(10)=IDATA1(12)
      IQTQ(11)=IDATA1(14)
      GOTO 800
C---- FILL IN *GENERAL CASE* AND IN IQTQ(13) A BLANK
  750 DO 760 I=1,3
  760 IQTQ(I+9)=IDATA1(I+23)
      IQTQ(13)=IDATA1(27)
C---- FILL IN *GAUSSIAN DISTRIBUTION*
      DO 770 I=1,2
      IQTQ(I+13)=IDATA1(I+14)
      IQTQ(I+15)=IDATA1(I+19)
  770 IQTQ(I+17)=IDATA1(I+21)
C---- CHANGE TO LORENTZIAN IF NECESSARY
      IF(NDISTR)800,800,780
  780 DO 790 I=1,7
  790 IQTQ(I+12)=IDATA1(I+16)
  800 CONTINUE
C---- OKAY, PRINT NOW EXTINCTION TYP SELECTED.
      WRITE(IOUT,20)(IQTQ(I),I=1,19)
C---- NOW GO ON WITH PARAMETERS, IF EXTINCTION IS COMPUTED.
      IF(ITYP.EQ.1)GOTO 1010
C---- PRINT HEADER
      WRITE(IOUT,22)
C---- FILL IN *EXTINCTION PARAMETERS X 10**-4 *
      DO 810 I=1,2
      IQTQ(I)=IDATA1(I+1)
      IQTQ(I+2)=IDATA1(I+27)
      IQTQ(I+4)=IDATA1(I+29)
  810 IQTQ(I+6)=IDATA1(I+31)
C---- PRINT THIS AND PARAMETERS.
      GOTO(1010,820,830,840,850,850,860),ITYP
  820 WRITE(IOUT,24)(IQTQ(I),I=1,8),BDF(3)
      GOTO 870
  830 WRITE(IOUT,26)(IQTQ(I),I=1,8),BDF(4)
      GOTO 870
  840 WRITE(IOUT,28)(IQTQ(I),I=1,8),BDF(4),BDF(3)
      GOTO 870
  850 JUNK3=1
      JUNK4=6
      IF(ITYP.EQ.6) THEN
        JUNK3=7
        JUNK4=12
      END IF
      WRITE(IOUT,30)(IQTQ(I),I=1,8),(BDF(I+4),I=JUNK3,JUNK4)
      GOTO 870
  860 CONTINUE
      WRITE(IOUT,32)(IQTQ(I),I=1,8),BDF(4),(BDF(4+I),I=1,6)
  870 CONTINUE
C
C
C----------------------------------------------------------------------
      ENTRY EXTRA
C----------------------------------------------------------------------
C
C
C---- NOW RELOAD VALUES AND  COMPUTE DOMAINSIZES AND/OR  MOSAIC SPREAD.
C
      ITYP=IBDF(1)+1
      IDISTR=IBDF(2)
      JUNK7= 0
      GOTO(1010,890,890,890,880,880,880),ITYP
  880 CONTINUE
      JUNK5=1
      JUNK4=41
      IF(ITYP.EQ.6) JUNK4 = 47
      GOTO 540
  890 CONTINUE
      QTQ(55)=5.8186
      IF(IDISTR.NE.0)QTQ(55)=3.2828
      JUNK1=2
      IF(ITYP.EQ.4) JUNK1=4
      IF(ITYP.EQ.7) JUNK1=4
      GOTO(1010,900,900,900,910,910,910),ITYP
  900 JUNK3=40
      JUNK4=40
      IF(ITYP.NE.2) JUNK4=41
      IF(ITYP.EQ.3) JUNK3=41
      GOTO 920
  910 JUNK3=42
      JUNK4=44
      IF(ITYP.EQ.6) THEN
        JUNK3 = 48
        JUNK4 = 50
      END IF
  920 JUNK2=33
      IF(ITYP.EQ.3) JUNK2 = 41
      IF(ITYP.EQ.6) JUNK2 = 41
      DO 930 I=1,8
      IQTQ(I+59)=IDATA1(I+JUNK2)
      IQTQ(I+67)=IDATA1(I+41)
  930 CONTINUE
      GOTO(1010,940,940,940,970,970,970),ITYP
  940 DO 960 I=JUNK3,JUNK4
      IF(QTQ(I).GT.1.E-10)GOTO 950
      JUNK7=1
  950 IF(JUNK7.NE.1) QTQ(I-JUNK3+80) = QTQ(55)/QTQ(I)
      IF(ITYP.EQ.3)  QTQ(I-JUNK3+80) = 0.0001*QTQ(I)
      IF(I.NE.JUNK3) QTQ(I-JUNK3+80) = 0.0001*QTQ(I)
  960 CONTINUE
      IF(ITYP.EQ.4)WRITE(IOUT,34)(IQTQ(I+67),I=1,8),QTQ(81)
      IF(ITYP.NE.3.AND.JUNK7.EQ.0)WRITE(IOUT,34)
     & (IQTQ(I+59),I=1,8),QTQ(80)
      IF(ITYP.EQ.3)WRITE(IOUT,34)(IQTQ(I+59),I=1,8),QTQ(80)
      GOTO 1010
  970 WRITE(IOUT,42)QTQ(1),QTQ(5),QTQ(9)
      IF(QTQ(1).LT.0. .OR. QTQ(5).LT.0. .OR. QTQ(9).LT.0.) GOTO 1010
      IF(ITYP.EQ.6)  THEN
C---- DOMAINSIZE
        QTQ(1)=0.0001/SQRT(QTQ(1))
        QTQ(2)=0.0001/SQRT(QTQ(5))
        QTQ(3)=0.0001/SQRT(QTQ(9))
      ELSE
C---- MOSAIC SPREAD
        QTQ(1)=QTQ(55)*SQRT(QTQ(1))
        QTQ(2)=QTQ(55)*SQRT(QTQ(5))
        QTQ(3)=QTQ(55)*SQRT(QTQ(9))
      END IF
 1000 CONTINUE
      WRITE(IOUT,44)(IQTQ(I+59),I=1,8),(QTQ(I),I=1,3)
 1010 CONTINUE
      RETURN
      END
                      SUBROUTINE JACOBI(QTQ)
C
C     CALCULATION OF EIGENVALUES OF THE ORTHOGONAL EXTINCTION TENSOR.
C          PROGRAMME ADAPTED FROM XRAY SYSTEM  BY U. ZUCKER
C
      DIMENSION QTQ(100)
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'E.INC'
C
   10 FORMAT(43H0CAUTION - TROUBLES IN DIAGONALIZING TENSOR)
C
  300 QTQ(19)=0.0
      I=0
      DO 330 J=1,3
      DO 330 K=1,3
      I=I+1
      IF(J-K)310,320,310
  310 QTQ(I+9)=0.0
      QTQ(19)=QTQ(19)+QTQ(I)**2
      GOTO 330
  320 QTQ(I+9)=1.0
      QTQ(19)=QTQ(19)+QTQ(I)
  330 CONTINUE
      QTQ(19)=QTQ(19)/3.
      QTQ(20)=0.0001*QTQ(19)
  340 QTQ(19)=0.333333*QTQ(19)
      I10=0
  350 IT=0
      I10=I10+1
      IF(I10-10)370,370,360
  360 WRITE(IOUT,10)
      GOTO 500
  370 CONTINUE
      JUNK1=1
  380 JUNK1=JUNK1+1
      JUNK2=0
  390 JUNK2=JUNK2+1
      I=3*(JUNK1-1)+JUNK2
      IF(ABS(QTQ(I))-QTQ(19))470,400,400
  400 J=3*(JUNK1-1)+JUNK1
      K=3*(JUNK2-1)+JUNK2
      L=3*(JUNK2-1)+JUNK1
      IT=IT+1
      QTQ(21)=0.5*(QTQ(K)-QTQ(J))
      TEMP=SQRT(QTQ(L)**2+QTQ(21)**2)
      IF(QTQ(21))420,410,410
  410 QTQ(22)=-QTQ(L)/TEMP
      GOTO 430
  420 QTQ(22)=QTQ(L)/TEMP
  430 QTQ(22)=SIGN(AMIN1(ABS(QTQ(22)),1.0),QTQ(22))
      SINR=QTQ(22)/(SQRT(2.0*(1.0+SQRT(1.0-QTQ(22)**2))))
      SINSQ=SINR**2
      COSSQ=1.0-SINSQ
      COSR=SQRT(COSSQ)
      SINCOS=SINR*COSR
      QTQ(23)=QTQ(J)
      QTQ(24)=QTQ(K)
      QTQ(25)=QTQ(L)
      DO 460 M=1,3
      I1=3*(JUNK2-1)+M
      I2=3*(JUNK1-1)+M
      I3=3*(M-1)+JUNK2
      I4=3*(M-1)+JUNK1
      QTQ(26)=QTQ(I1+9)
      QTQ(I1+9)=QTQ(I1+9)*COSR-QTQ(I2+9)*SINR
      QTQ(I2+9)=QTQ(26)*SINR+QTQ(I2+9)*COSR
      IF(M-JUNK2) 440,460,440
  440 IF(M-JUNK1) 450,460,450
  450 QTQ(26)=QTQ(I3)
      QTQ(I3)=QTQ(I3)*COSR-QTQ(I4)*SINR
      QTQ(I4)=QTQ(26)*SINR+QTQ(I4)*COSR
      QTQ(I1)=QTQ(I3)
      QTQ(I2)=QTQ(I4)
  460 CONTINUE
      QTQ(K)=QTQ(24)*COSSQ+QTQ(23)*SINSQ-2.0*QTQ(25)*SINCOS
      QTQ(J)=QTQ(24)*SINSQ+QTQ(23)*COSSQ+2.0*QTQ(25)*SINCOS
      QTQ(L)=(QTQ(24)-QTQ(23))*SINCOS+QTQ(25)*(COSSQ-SINSQ)
      QTQ(I)=QTQ(L)
  470 IF(JUNK2-JUNK1+1)390,480,480
  480 IF(JUNK1-3)380,490,490
  490 IF(IT)500,500,350
  500 IF(QTQ(19)-QTQ(20))510,510,340
  510 RETURN
C
      END
