                           SUBROUTINE INPUT (IEF,                               
     &INEUT,LPCORR,IFOBS,HELPE,HELPF,COS2TM,IEQPMT,ISKI,IVW,JJJ,IPSI)     
C                                                                               
C     SUBROUTINE TO READ STRUCTURE FACTOR CARDS, COMPUTE SCATTERING FAC-        
C     TORS , POLARISATION  CORRECTION,  LORENTZCORRECTION   AND WEIGHTS.        
C     ROUTINE IS IN CASE OF ANISOTROPIC EXTINCTION ROUTINE CREC CALLING,        
C                      WHICH COMPUTES VECTORS USED.                             
C                                                                               
C                          WRITTEN BY UDO ZUCKER                                
C               W.F. KUHS GOT VECTOR BUSINESS RIGHT May 1988
C               W.F. KUHS ADDED OPTION FOR PSI.NE.0 Feb 1992
C                                                                               
C     INPUT VARIABLES                                                           
C     ---------------                                                           
C                                                                               
C     IEF      =  1,2,3,4,5 = MOD(NEW(8),5)     GIVES  TYPE OF HKL-CARD.        
C     INEUT    =  0/1                           XRAY/NEUTRON  RADIATION.        
C     LPCORR   =  0/1  DO NOT/ DO  A LORENTZ (POLARISATION)  CORRECTION.        
C     IFOBS    =  0/1  INTENSITY/SCATTERING FACTORS        ON  HKL-CARD.        
C     HELPE    =  COEEFICIENT OF POLARISATION CORRECTION,  COMPUTED  IN         
C                 PROGRAMME REFINE.                                             
C     HELPF    =  LIKE HELPE                                                    
C     COS2TM   =  COS ( ANGLE OF MONOCHROMATOR*2 )**2 
C     IPSI     =  0/1  NO PSI-VALUES / PSI-VALUES ON HKL-CARDS
C                                                                               
      DIMENSION HK(3),SAVE(60),KDEGN(2640),AA(2640),ISH(3,120),BBETA(6,6        
     &0)                                                                         
C                                                                               
      COMMON/D/ACALC,AIND(60),ANOM(2,20),BCALC,BIND(60),DADB(6,60),DADO(        
     & 60),DADX(3,60),DBDB(6,60),DBDO(60),DBDX(3,60),H(3),IH(3),ISCAT(2,        
     & 60), ISCL, MFSEL,OCCA(60), OCCUP(60),SCALE(10),SCAT(20),SITE(60),        
     & DADC(10,60), DBDC(10,60), DADD(15,60), DBDD(15,60),AMINUS,BMINUS,        
     & DRDE(7),NDEGN,VORFAK,IPAREX(7),DAALFA(60),DBALFA(60)                     
      COMMON/TAPE/IN,INN,IOUT,ISC1,ISC2,ISC3,ISC4,ISC5,IPUN                     
      COMMON/E/DMAX,DMIN,RECELL(6),EXANG(6),EXTRAN(12)                          
      COMMON/F/DELTA,FCALC,FOBS,FOMAX,FOMIN,IREJ,MREF,MREJ,NREJ,NUMBER(3        
     & 6), RDEN(36), RHO,RHOMN, RHOMX,RNUM(36),WT,EXBETA,EXX,EXY,GAMMAP,        
     & FISGN,IFSGN,EXTIN(6),KEXTRF(7),DELMAX,IANIGC,GAMMAS,IZACH,RNEUEX(        
     & 7),NEUEXT,NDISTR,NIEXT,WTA(3),WTB(3),RHODIV                              
      COMMON/G/BUFF(1023),FCURVE(9,20),IFMT(18),IWT,NSCAT,BETFAC,WVLNTH         
     &,WTFAC                                                                    
      COMMON/POL/INTER(20),X(35,4),Y(35,4),INTPOL                               
      COMMON/H/UB(9),CHI,PHI,TH,UBINV(9),TBAR,DIN(3),DOUT(3),IEQM,BBETA         
      COMMON/FOUR/WINKEL(3),IUB,AG,BG,CG                                        
C                                                                               
      EQUIVALENCE (SAVE,AIND),(KDEGN,AA,DADC(1,37)),(ISH,BBETA)                 
C                                                                               
 5000 FORMAT(/13H0ERRORMESSAGE/1X,12(1H-)/1H )                                  
 5001 FORMAT(55H0AN AVERAGED MEAN PATH LENGTH LESS THAN 1.E-8 CMS ENCOU,        
     & 16HNTERED ON INPUT./30H0THAT SEEMS TO BE MEANINGLESS./                   
     & 61H0CHECK THE CONSISTENCY OF THE PARAMETER NEW(8) AND THE HKL-CA,        
     & 12HRDS READ IN./16H0JOB TERMINATED.)                                     
 5002 FORMAT(55H0A NULL SCALINGGROUP INDICATOR ENCOUNTERED ON INPUT.   /        
     & 61H0CHECK THE CONSISTENCY OF THE PARAMETER NEW(8) AND THE HKL-CA,        
     & 12HRDS READ IN./16H0JOB TERMINATED.)                                     
 5003 FORMAT(55H0A NULL STANDARD DEVIATION OR WEIGHT ENCOUNTERED ON INP,        
     & 3HUT./30H0THAT SEEMS TO BE MEANINGLESS./16H0JOB TERMINATED.)             
 5004 FORMAT(55H0ALL HKL INDICES OF A REFLECTION READ IN ARE ZERO.     /        
     & 30H0THAT SEEMS TO BE MEANINGLESS./16H0JOB TERMINATED.)                   
 5005 FORMAT(55H0A REJECTION FLAG INDICATOR WITH ILLEGAL VALUE ENCOUNTE,        
     & 16HRED ON INPUT    /38H0ONLY VALUES OF ONE AND TWO ARE LEGAL./           
     & 16H0JOB TERMINATED.)                                                     
C                                                                               
C---- DECLARE SYMBOL FOR END OF HKL CARDS.  (END MUST BE GIVEN BY THE           
C     FIRST THREE CHARACTERS OF A HKL CARD)                                     
C                                                                               
      DATA IEND/3HEND/                                                          
C                                                                               
C---- STORE   BETFAC(TBAR=CONSTANT) AND GAMMAS(TBAR=CONSTANT) TO KNOW           
C     THESE VALUES FOR UNFIXED TIBAR.                                           
C                                                                               
      BETFA=BETFAC                                                              
      GAMMS=GAMMAS                                                              
C                                                                               
C---- BRANCH IN ACCORDANCE TO IEF=MOD(NEW(8),5), FIRST CHECK ON PSI
C
      IF(IPSI.EQ.1) GOTO 6 
C                                                                               
      GO TO (1,2,3,4,5),IEF                                                     
C                                                                               
C---- THERE IS AN AVERAGED TIBAR FOR ALL REFLECTIONS ( READ IN REFINE ).        
C                                                                               
    1 READ(INN,IFMT,END=120)LABEL,IH,MREJ,FOBS,SIGMA                     
C                                                                               
C---- END OF LIST READ?                                                         
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
C                                                                               
C---- I IS SCALEFACTOR NUMBER FIXED TO I+1 = 1                                  
      I=0                                                                       
      GO TO 1000                                                                
C                                                                               
C---- TIBAR IN CENTIMETERS ON EACH HKL-CARD FOR EACH REFLECTION.                
C                                                                               
    2 READ(INN,IFMT,END=120)LABEL,IH,MREJ,FOBS,TBAR,SIGMA                  
C                                                                               
C---- BETFAC, GAMMAS  ARE  COEFFICIENTS IN EXTINCTION CORRECTION.               
C     THESE HAVE TO BE MULTIPLIED BY TIBAR, HERE WRITTEN AS TBAR.               
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
      IF(TBAR.LT.1.E-8) GOTO 6001                                               
      BETFA=BETFAC*TBAR                                                         
      GAMMS=GAMMAS*TBAR                                                         
      I=0                                                                       
      GOTO 1000                                                                 
C                                                                               
C---- AVERAGED TIBAR AND SCALEGROUP INDICATOR ON CARD.                          
C                                                                               
    3 READ(INN,IFMT,END=120)LABEL,IH,MREJ,I,FOBS,SIGMA                  
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
      IF(I.EQ.0) GOTO 6002                                                      
      GOTO 1000                                                                 
C                                                                               
C---- FOR EACH REFLECTION TIBAR AND SCALEGROUP ON CARD.                         
C                                                                               
    4 READ(INN,IFMT,END=120)LABEL,IH,MREJ,I,FOBS,TBAR,SIGMA                 
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
      IF(TBAR.LT.1.E-8) GOTO 6001                                               
      IF(I.EQ.0) GOTO 6002                                                      
      BETFA=BETFAC*TBAR                                                         
      GAMMS=GAMMAS*TBAR                                                         
      GOTO 1000                                                                 
C                                                                               
C---- ONLY HKL READ FOR STRUCTUREFACTOR CALCULATION ONLY.                       
    5 READ(INN,IFMT,END=120)LABEL,IH                                            
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
C                                                                               
C---- SET REST OF INFO. TO SOME VALUES.                                         
      FOBS = 0.0                                                                
      MREJ = 1                                                                  
      I    = 0                                                                  
      SIGMA= 1.0                                                                
      GOTO 1000
C
C---- HERE WE GO FOR READING PSI-VALUES
    6 GO TO (31,32,33,34,35),IEF                                              
C                                                                               
C---- THERE IS AN AVERAGED TIBAR FOR ALL REFLECTIONS ( READ IN REFINE ).        
C                                                                               
   31 READ(INN,IFMT,END=120)LABEL,IH,MREJ,FOBS,SIGMA,PSI                      
C                                                                               
C---- END OF LIST READ?                                                         
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
C                                                                               
C---- I IS SCALEFACTOR NUMBER FIXED TO I+1 = 1                                  
      I=0                                                                       
      GO TO 1000                                                                
C                                                                               
C---- TIBAR IN CENTIMETERS ON EACH HKL-CARD FOR EACH REFLECTION.                
C                                                                               
   32 READ(INN,IFMT,END=120)LABEL,IH,MREJ,FOBS,TBAR,SIGMA,PSI                  
C                                                                               
C---- BETFAC, GAMMAS  ARE  COEFFICIENTS IN EXTINCTION CORRECTION.               
C     THESE HAVE TO BE MULTIPLIED BY TIBAR, HERE WRITTEN AS TBAR.               
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
      IF(TBAR.LT.1.E-8) GOTO 6001                                               
      BETFA=BETFAC*TBAR                                                         
      GAMMS=GAMMAS*TBAR                                                         
      I=0                                                                       
      GOTO 1000                                                                 
C                                                                               
C---- AVERAGED TIBAR AND SCALEGROUP INDICATOR ON CARD.                          
C                                                                               
   33 READ(INN,IFMT,END=120)LABEL,IH,MREJ,I,FOBS,SIGMA,PSI                   
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
      IF(I.EQ.0) GOTO 6002                                                      
      GOTO 1000                                                                 
C                                                                               
C---- FOR EACH REFLECTION TIBAR AND SCALEGROUP ON CARD.                         
C                                                                               
   34 READ(INN,IFMT,END=120)LABEL,IH,MREJ,I,FOBS,TBAR,SIGMA,PSI                 
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
      IF(TBAR.LT.1.E-8) GOTO 6001                                               
      IF(I.EQ.0) GOTO 6002                                                      
      BETFA=BETFAC*TBAR                                                         
      GAMMS=GAMMAS*TBAR                                                         
      GOTO 1000                                                                 
C                                                                               
C---- ONLY HKL READ FOR STRUCTUREFACTOR CALCULATION ONLY.                       
   35 READ(INN,IFMT,END=120)LABEL,IH                                            
C                                                                               
      IF(LABEL.EQ.IEND) GOTO 120                                                
C                                                                               
C---- SET REST OF INFO. TO SOME VALUES.                                         
      FOBS = 0.0                                                                
      MREJ = 1                                                                  
      I    = 0                                                                  
      SIGMA= 1.0                                                                
C                                                                               
C---- SKIP UNWANTED REFLECTIONS READ IN REFIN FROM UNIT 16                      
C                                                                               
 1000 IF(ISKI.NE.1) GOTO 1001                                                   
      DO 804 III=1,JJJ                                                          
      DO 803 K=1,3                                                              
      IF(IH(K).NE.ISH(K,III)) GOTO 804                                          
  803 CONTINUE                                                                  
      IVW = 1                                                                   
      RETURN                                                                    
  804 CONTINUE                                                                  
C                                                                               
C---- REFLECTION PROCESSING STARTS.       (RHO=SIN(THETA)/LAMBDA**2)            
C                                                                               
C     IN ACCORDANCE TO MODEL OF EXTINCTION COEFFICIENTS ARE SELECTED.           
C                                                                               
 1001 IF(IZACH.EQ.0)EXBETA=GAMMS                                                
C                                                                               
C---- NEGATIVE IINTENSITIES OR STRUCTUREFACTORS ARE ILLEGAL.                    
C                                                                               
      IF(FOBS.LT.0.0)FOBS=0.0                                                   
C                                                                               
C---- CALCULATE HERE REJECTION WITH RESPECT TO SIGMA                            
C                                                                               
      IF(IREJ.NE.2.OR.DELMAX.EQ.0.0) GOTO 10                                    
      IF(FOBS/SIGMA.LT.DELMAX) MREJ = 2                                         
      IF(FOBS/SIGMA.GE.DELMAX) MREJ = 1                                         
C                                                                               
C---- INTENSITY OR AMPLITUDE ON CARD?                                           
C                                                                               
   10 IF(IFOBS.NE.0) GOTO 11                                                    
      FOBS=SQRT(FOBS)                                                           
   11 RHO=0.0                                                                   
C                                                                               
C---- LTHAN DESIGNATOR OF HKL-CARDS AND SO OF THE XRAY SYSTEM MUST DE-          
C     CREASED BY ONE FOR THIS PROGRAMME.                                        
C                                                                               
      IF(MREJ.EQ.0) MREJ=1                                                      
      MREJ=MREJ-1                                                               
C                                                                               
C---- ALL MILLER INDICES CANNOT BE ZERO.                                        
C                                                                               
      IEF=ABS(IH(1))+ABS(IH(2))+ABS(IH(3))                                      
      IF(IEF)13,12,13                                                           
   12 GOTO 6004                                                                 
C                                                                               
C---- REJECTION FLAG HAVE TO BE ONE ORE TWO.                                    
C                                                                               
   13 IF(MREJ.EQ.0 .OR. MREJ.EQ.1) GOTO 15                                      
      GOTO 6005                                                                 
C                                                                               
C---- SIGMA MUST GREATER ZERO. (PROBLEM OF WEIGHTING SCHEME.)                   
   15 IF(SIGMA)16,14,16                                                         
   14 GOTO 6003                                                                 
                                                                                
C                                                                               
C---- CALCULATE (SIN(THETA)/LAMBDA)**2                                          
C                                                                               
   16 DO 20 K=1,3                                                               
      H(K)=FLOAT(IH(K))                                                         
      HK(K)=H(K)*RECELL(K)                                                      
   20 RHO=RHO+HK(K)**2                                                          
      RHO=(RHO+2.0*(HK(1)*HK(2)*RECELL(6)+HK(1)*HK(3)*RECELL(5)+HK(2)*HK        
     & (3)*RECELL(4)))*0.25                                                     
      SIG=SQRT(RHO)                                                             
C                                                                               
C---- LORENTZ/(POLARISATIONS) COORECTION OF GEOMETRY OF AN EULER CRADLE.        
C     IN CASE OF NEUTRON, ONLY LORENTZ CORRECTION WILL BE DONE.                 
C                                                                               
      IF(LPCORR.EQ.0) GOTO 2000                                                 
      SINT2=SIN(2.*ASIN(SIG*WVLNTH))                                            
      IF(INEUT) 1997,1997,1998                                                  
 1997 COS2T2=1.-SINT2**2                                                        
      RLP=(1.+COS2TM)/(HELPF+HELPE*COS2T2)*SINT2                                
      GOTO 1999                                                                 
 1998 RLP=SINT2                                                                 
 1999 SIGMA=SIGMA*RLP                                                           
      RLP=SQRT(RLP)                                                             
      FOBS=FOBS*RLP                                                             
C                                                                               
C---- SET FOBS=0 TO A SMALL VALUE.                                              
 2000 IF(FOBS.GT.0.00001)GOTO 2001                                              
      FOBS=1.E-7                                                                
C                                                                               
C---- COMPUTE MINIMUM AND MAXIMUM OF FOBS AND OF  SIN(T)/LAMBDA.                
C                                                                               
 2001 IF(FOBS-FOMAX)22,22,21                                                    
   21 FOMAX=FOBS                                                                
   22 IF(FOBS-FOMIN)23,24,24                                                    
   23 FOMIN=FOBS                                                                
   24 IF(RHO-RHOMX)26,26,25                                                     
   25 RHOMX=RHO                                                                 
   26 IF(RHO-RHOMN)27,28,28                                                     
   27 RHOMN=RHO                                                                 
C                                                                               
C---- CALCULATE EFFECT OF POLARISATION AND SCATT.-ANGLE ON EXTINCTION.          
C                                                                               
   28 IF (IZACH.NE.1) GOTO 1010                                                 
C                                                                               
C---- FOR ZACHARIASEN CORRECTION.                                               
C                                                                               
C---- ANGLE.                                                                    
      THETA2=2*ASIN(SIG*WVLNTH)                                                 
      EXBETA=BETFA/SIN(THETA2)                                                  
      IF(INEUT.NE.0) GOTO 1010                                                  
C---- POLARISATIONCORRECTION.   (XRAY ONLY)                                     
      COS2T2=COS(THETA2)**2                                                     
      Q=HELPE*COS2T2                                                            
      Q=(HELPF+Q)/(HELPF+Q*COS2T2)                                              
      EXBETA=EXBETA/Q                                                           
C                                                                               
C---- ANISOTROPIC EXTINCTION? NO  GO TO 29                                      
C                                                                               
 1010 IF(IUB.NE.1) GOTO 29                                                      
C---- CALCULATE FROM  UB-MATRIX AND HKL-VALUE SETTING ANGLES OF DIFFR.          
C                                                                               
      CALL CREC(IZACH,WVLNTH,IEQPMT,PSI)                                      
C                                                                               
      IF(IZACH.NE.1)GO TO 1011                                                  
C                                                                               
C---- SETTING ANGLE OF DIFFRACTOMETER FOR EXTINCTION OF COPPENS AND             
C     & HAMILTON .(FROM NOW BUSING AND LEVY GEOMETRY)                           
C                                                                               
      EXANG(1)=-CHI                                                             
      EXANG(2)=PHI                                                              
      EXANG(3)=TH                                                               
      GO TO 29                                                                  
C                                                                               
C---- BECKER & COPPENS. UNITVECTORS IN DIRECTION OF DIFFRACTED BEAM .           
C     (REAL LATTICE AS BASE)                                                    
C                                                                               
 1011 EXANG(1)=DOUT(1)
      EXANG(2)=DOUT(2)
      EXANG(3)=DOUT(3)
C                                                                               
C---- UNITVECTOR NORMAL TO THAT PLANE, IN WHICH INCIDENT AND                    
C     REFLECTED BEAM LAYING. (BASE IS REAL CELL)                                
C                                                                               
      EXANG(4)=DIN(1)
      EXANG(5)=DIN(2)
      EXANG(6)=DIN(3)
C                                                                               
   29 DO 40 L=1,NSCAT                                                           
C---- COMPUTE VALUE OF SCATTERING FACTORS OF THIS REFECTION.                    
C                                                                               
      IF(INTER(L))41,42,41                                                      
C---- GIVEN IN TABULATION FORM.                                                 
   41 LL=L                                                                      
      CALL INPOL(RHO,LL,RESULT)                                                 
      SCAT(L)=RESULT+ANOM(1,L)                                                  
      GO TO 40                                                                  
C---- GIVEN IN EXPONENTIAL FORM.                                                
   42 SCAT(L)=FCURVE(9,L)+ANOM(1,L)                                             
      DO 43 M=1,8,2                                                             
      IF(FCURVE(M+1,L)*RHO-88..GT.0.)GOTO 43                                    
      SCAT(L)=SCAT(L)+FCURVE(M,L)*EXP(-FCURVE(M+1,L)*RHO)                       
   43 CONTINUE                                                                  
   40 CONTINUE                                                                  
C                                                                               
C---- COMPUTE WEIGHT.                                                           
C                                                                               
      WT=SIGMA                                                                  
C                                                                               
      GO TO (50,90,53,57,58,59,45,55,53,45),IWT                                 
C                                                                               
C---- IWT = 1 SIGMA(F*F) ON CARD, COMPUTE SIGMA(F), USE FORMULA OF REES.        
C         = 2 SIGMA(F)   ON CARD.                                               
C         = 3 WEIGHT     ON CARD. IF REFINEMENT ON F*F SELECTED,  WEIGHT        
C                                 OF F*F, OTHERWISE FOR F.                      
C         = 4 LIKE 1, BUT WITHOUT FORMALISM OF REES.                            
C         = 5 SIGMA(F*F) ON CARD USE ROUTINE WEIGHT.                            
C         = 6 SIGMA(F) ON CARD USE SPECIAL ROUTINE WEIGHT.                      
C         = 7 LIKE 1, BUT ADD AN IGNORANCE FACTOR.                              
C         = 8 LIKE 2, BUT ADD AN IGNORANCE FACTOR.                              
C         = 9 LIKE 3, BUT ADD AN IGNORANCE FACTOR.                              
C         =10 LIKE 4, BUT ADD AN IGNORANCE FACTOR.                              
C                                                                               
C---- USE FORMULA BY REES FOR CALCULATION OF SIGMA(F) (ACTA CRYST.              
C     A32, 483-488 (1976)                                                       
C                                                                               
   45 WT = SQRT(WT*WT + (FOBS**2 * WTFAC)**2 )                                  
      IF(IWT.EQ.10) GOTO 57                                                     
C                                                                               
   50 IF(MREF) 502,501,502
  501 IF(FOBS*FOBS-WT) 52,52,51                                                 
   51 WT=0.5*WT/FOBS                                                            
      GO TO 90                                                                  
   52 WT=SQRT(WT)*0.5                                                           
      GOTO 90                                                                   
  502 WT = 0.5*WT/FOBS
      GOTO 90
C                                                                               
   53 WT=1.0/SQRT(WT)                                                           
      IF(MREF)54,56,54                                                          
   54 IF(IWT.EQ.6) GOTO 45                                                      
C---- WT IS HERE SIGMA(I)                                                       
      WT=WT*0.5/FOBS                                                            
      GO TO 90                                                                  
C                                                                               
   55 WT=SQRT(WT*WT+(WTFAC*FOBS)**2)                                            
      GOTO 90                                                                   
C                                                                               
   56 IF(IWT-6)90,55,90                                                         
C                                                                               
   57 WT=0.5*SIGMA/FOBS                                                         
      GOTO 90                                                                   
C                                                                               
C  58 WT = 0.5*SIGMA/FOBS                                                       
   58 CALL WEIGHT                                                               
      GOTO 90                                                                   
C                                                                               
   59 CALL WEIGHT                                                               
C                                                                               
C---- SET SCALE FACTOR FOR REFLECTION                                           
C                                                                               
   90 IF(I)100,110,100                                                          
  100 ISCL=I                                                                    
      GOTO 120                                                                  
C                                                                               
C---- ONLY ONE SCALE GROUP                                                      
  110 ISCL=1                                                                    
      GOTO 120                                                                  
C                                                                               
C---- NO MORE REFLECTION                                                        
C                                                                               
  120 IEF=0                                                                     
      IF(LABEL.EQ.IEND) IEF=1                                                   
C                                                                               
      RETURN                                                                    
C                                                                               
C---- OUTPUT PART FOR ERRORS DETECTED ON HKL-CARDS.                             
C                                                                               
 6001 WRITE(IOUT,5000)                                                          
      WRITE(IOUT,5001)                                                          
      CALL EXIT                                                                 
 6002 WRITE(IOUT,5000)                                                          
      WRITE(IOUT,5002)                                                          
      CALL EXIT                                                                 
 6003 WRITE(IOUT,5000)                                                          
      WRITE(IOUT,5003)                                                          
      CALL EXIT                                                                 
 6004 WRITE(IOUT,5000)                                                          
      WRITE(IOUT,5004)                                                          
      CALL EXIT                                                                 
 6005 WRITE(IOUT,5000)                                                          
      WRITE(IOUT,5005)                                                          
      CALL EXIT                                                                 
C                                                                               
      RETURN                                                                    
      END                                                                       
                 SUBROUTINE CREC(IZACH,WVLNTH,IEQPMT,PSI)                       
C                                                                               
C     COMPUTATION OF DIFFRACTOMETER SETTING ANGLES FROM  UB-MATRIX  AND         
C     FROM HKL-VALUE AS WELL AS CALCULATION OF SOME  UNIT VECTORS  USED         
C                  IN ANISOTROPIC   EXTINCTION THEORY.                          
C                                                                               
C                                                                               
C     WRITTEN BY UDO ZUCKER,USING THE PHILIPS PW1100 PROGRAMME AS HELP.         
C   AND CORRECTED TO GET VECTORS RIGHT FOR BUSING & LEVY BY W.F.KUHS OCT 84     
C     PSI.NE.0  and  PICKER-STOE added by W.F.KUHS FEB 92               
C                       
C     INPUT VARIABLES                                                           
C     ---------------                                                           
C                                                                               
C     IZACH = 0/1         BECKER&COPPENS /ZACHARIASEN EXTINCTION THEORY.        
C     WVLNTH= WAVELENGTH                                                        
C     IEQPMT= FLAG FOR GEOMETRY OF DIFFRACTOMETER. (IN THE CURRENT  VER-        
C             SION IS   1 = PW1100 (PHILIPS), 2 = BUSING & LEVY GEOMETRY        
C                       3 = PICKER-STOE, 4 = SYNTEX P3
C---- BUSING & LEVY GEOMETRY CORRECTED 27.01.84 W.F.KUHS                        
C                                                                               
      DIMENSION RECCO(6)                                                        
C                                                                               
      COMMON/H/UB(9),CHI,PHI,TH,UBINV(9),TBAR,DIN(3),DOUT(3)                    
      COMMON/D/ACALC,AIND(60),ANOM(2,20),BCALC,BIND(60),DADB(6,60),DADO(        
     & 60),DADX(3,60),DBDB(6,60),DBDO(60),DBDX(3,60),H(3),IH(3),ISCAT(2,        
     & 60),ISCL,MFSEL,OCCA(60),OCCUP(60), SCALE(10), SCAT(20), SITE(60),        
     & DADC(10,60),DBDC(10,60), DADD(15,60), DBDD(15,60), AMINUS,BMINUS,        
     & DRDE(7),NDEGN,VORFAK,IPAREX(7),DAALFA(60),DBALFA(60)                     
C                                                                               
      GOTO (100,100,300,400) , IEQPMT
C
C---- PHILIPS DIFFRACTOMETER PW1100
C
  100 DO 10 I=1,3                                                               
      RECCO(I)=0                                                                
        DO 10 J=1,3                                                             
        RECCO(I)=RECCO(I)+UB(3*(I-1)+J)*IH(J)                                   
   10 CONTINUE                                                                  
C                                                                               
      T1=RECCO(1)*RECCO(1)+RECCO(2)*RECCO(2)                                    
      T2=T1+RECCO(3)*RECCO(3)                                                   
C---- TH IS THETA IN RADIAN.                                                    
      TH=ASIN(0.5*SQRT(T2))                                                     
C---- PHI,CHI IS CLEAR, IS NOT IT.                                              
      IF(RECCO(1))1040,1030,1040                                                
 1030 IF(RECCO(2))1040,1035,1040                                                
 1035 PHI=1.5707963                                                     
      GO TO 1045                                                                
 1040 PHI=ATAN2(RECCO(1),RECCO(2))                                              
 1045 CHI=-ATAN2(RECCO(3),SQRT(T1))                                             
C                                                                               
C---- TRANSFORMATION OF DIFFRACTOMETER SYSTEMS.                                 
C                                                                               
      GOTO (1,2,3,3),IEQPMT                                                       
C                                                                               
C---- PHILIPS DIFFRACTOMETER PW1100                                             
C                                                                               
    1 GOTO 4                                                                    
C---- BUSING AND LEVY GEOMETRY                                                  
C                                                                               
    2 CHI = -CHI                                                               
      PHI = 1.5707963 - PHI                                                     
      GOTO 4                                                                    
C---- FREE FOR USER DIFFRACTOMETER                                              
    3 CALL EXIT                                                                 
C     GOTO 4                                                                    
C                                                                               
    4 IF (IZACH.EQ.1) RETURN                                                    
C                                                                               
C---- GO ON WITH BECKER AND COPPENS. HERE AS BASE THE REAL CRYSTAL CELL         
C     IS USED.  (COMPUTATION LIKE XRAY 80, PSI-ANGLE IS FIXED.)                 
C                                                                               
      SCHI = SIN(CHI)                                                           
      CCHI = COS(CHI)                                                           
      STH  = SIN(TH)                                                            
      CTH  = COS(TH)                                                            
      SOM =  STH
      COM =  CTH
      CTTH = COS(2.*TH)
      STTH = SIN(2.*TH)
      SPHI = SIN(PHI)                                                           
      CPHI = COS(PHI)                                                           
      IF (IEQPMT.EQ.2) GOTO 201
C                                                                               
C---- VECTOR IN DIRECTION OF INCIDENT BEAM (IN THE PHI SYSTEM)                  
C                                                                               
      RECCO(1) = -SPHI*CCHI*STH+CPHI*CTH                                        
      RECCO(2) = -CPHI*CCHI*STH-SPHI*CTH                                        
      RECCO(3) = SCHI*STH                                                       
C                                                                               
C---- VECTOR IN DIRECTION OF REFLECTED BEAM (PHI SYSTEM)                        
C                                                                               
      RECCO(4) = SPHI*CCHI*STH+CPHI*CTH                                         
      RECCO(5) = CPHI*CCHI*STH-SPHI*CTH                                         
      RECCO(6) = -RECCO(3)
      GOTO 202                                                      
C
C---- CALCULATE VECTOR OF BEAM IN AND NORMAL TO BEAM OUT IN SCATTERING
C     PLANE IN PHI-SYSTEM
C
  201 RECCO(1) =  SOM*CCHI*CPHI + COM*SPHI
      RECCO(2) = -SOM*CCHI*SPHI + COM*CPHI
      RECCO(3) = -SOM*SCHI
      RECCO(4) =  COM*CCHI*CPHI - SOM*SPHI
      RECCO(5) = -COM*CCHI*SPHI - SOM*CPHI
      RECCO(6) = -COM*SCHI
C                                                                               
C---- SWITCH ACCORDING TO DIFFRACTOMETER SYSTEM                                 
C                                                                               
  202 GOTO (6,7,8,8),IEQPMT                                                       
C                                                                               
    6 GOTO 9                                                                    
C                                                                               
C   7 RECCS = RECCO(1)                                                          
    7 RECCO(1) = -RECCO(1)                                                       
C     RECCO(2) = RECCO(2)                                                          
C     RECCO(3) = RECCO(3)                                                      
C                                                                               
C     RECCS = RECCO(4)                                                          
      RECCO(4) = -RECCO(4)                                                       
C     RECCO(5) = RECCO(5)                                                          
C     RECCO(6) = RECCO(6)                                                      
C     WRITE(15,9998) RECCO
C9998 FORMAT(' RECCO :',1X,6F9.4)
      GOTO 9                                                                    
C                                                                               
    8 CALL EXIT                                                                 
C                                                                               
C---- TRANSFORMATION IN BASE OF RECIPROCAL SYSTEM.                              
C                                                                               
    9 JUNK=0                                                                    
C                                                                               
      GOTO (203,204,8,8),IEQPMT  
C---- BUSING & LEVY CASE
  204 DO 1065 J=1,3                                                             
      DIN(J)=0.0                                                                
      DOUT(J)=0.0                                                               
      DO 1060 K=1,3                                                             
      JUNK   = JUNK+1                                                           
      DIN(J) = DIN(J)+UBINV(JUNK)*RECCO(K)*WVLNTH                               
      DOUT(J)= DOUT(J)+UBINV(JUNK)*RECCO(3+K)*WVLNTH                            
 1060 CONTINUE                                                                  
      DOUT(J) = STTH*DOUT(J) - CTTH*DIN(J)
 1065 CONTINUE
      GOTO 20
C---- PHILIPS CASE
  203 DO 1066 J=1,3                                                             
      DIN(J)=0.0                                                                
      DOUT(J)=0.0                                                               
      DO 1061 K=1,3                                                             
      JUNK   = JUNK+1                                                           
      DIN(J) = DIN(J)+UBINV(JUNK)*RECCO(K)*WVLNTH                               
      DOUT(J)= DOUT(J)+UBINV(JUNK)*RECCO(3+K)*WVLNTH                            
 1061 CONTINUE                                                                  
 1066 CONTINUE
C     WRITE(15,9999) IEQPMT,IH,TH,OM,CHI,PHI,DIN,DOUT                                 
C9999 FORMAT(1X,I2,4X,3I4,1X,4F8.4,2X,3F8.4,2X,3F8.4)                                 
      GOTO 20                                                                  
C
C---- PICKER-STOE DIFFRACTOMETER
C
C---- FIRST CALCULATE BISECTING ANGLES
C
  300 DO 310 I=1,3                                                             
      RECCO(I)=0                                                                
        DO 310 J=1,3                                                           
        RECCO(I)=RECCO(I)+UB(3*(I-1)+J)*IH(J)                                   
  310 CONTINUE                                                                
C                                                                               
      T1=RECCO(1)*RECCO(1)+RECCO(2)*RECCO(2)                                    
      T2=T1+RECCO(3)*RECCO(3)                                                   
C---- TH IS THETA IN RADIAN.                                                    
      TH=ASIN(0.5*SQRT(T2))                                                     
C---- PHI,CHI IS CLEAR, IS IT NOT                                          
      IF(RECCO(1))3040,3030,3040                                                
 3030 IF(RECCO(2))3040,3035,3040                                                
 3035 PHI=1.5707963                                                            
      GO TO 3045                                                                
 3040 PHI=-ATAN2(RECCO(2),RECCO(1))                                              
 3045 CHI=ATAN2(RECCO(3),SQRT(T1))
C     WRITE(15,9998)IH,TH,OM,CHI,PHI,PSI
C9998 FORMAT(1X,3I4,1X,5F8.4)                                             
C---- NOW CALCULATE ANGLES FOR PSI.NE.O
C---- FIRST STORE BISECTING VALUES
      PHIZ = PHI
      CHIZ = CHI
      PSI = PSI*0.017453      
      IF(PSI.EQ.0) GOTO 3046
C---- CALCULATE OMEGA FOR PSI-ROTATION
      IF(ABS(ABS(CHIZ)-1.5707963).LT.1.E-8) GOTO 3046
      OM = ATAN(SIN(PSI)/TAN(CHIZ))
      COSOM = COS(OM)
      GOTO 3047
 3046 COSOM = 1.0
      OM = 0.0
 3047 SINOM = SIN(OM)
C---- CALCULATE TRUE PHI AND CHI
      CHI = ATAN2(SIN(CHIZ)/COSOM,COS(CHIZ)*COS(PSI))
      DELPHI = ATAN2(-SINOM/COS(CHIZ),COSOM*COS(CHI)/COS(CHIZ))
      PHI = PHIZ - DELPHI
C---- PREPARE FOR VECTOR CALCULATION
      SCHI = SIN(CHI)
      CCHI = COS(CHI)
      SPHI = SIN(PHI)
      CPHI = COS(PHI)
      SOM = SIN(TH-OM)
      COM = COS(TH-OM)
      STTH = SIN(2.*TH)
      CTTH = COS(2.*TH)
C---- CALCULATE VECTOR OF BEAM IN AND BEAM OUT IN PHI-SYSTEM
      RECCO(1) =  SOM*CCHI*CPHI + COM*SPHI
      RECCO(2) = -SOM*CCHI*SPHI + COM*CPHI
      RECCO(3) =  SOM*SCHI
      RECCO(4) =  COM*CCHI*CPHI - SOM*SPHI
      RECCO(5) = -COM*CCHI*SPHI - SOM*CPHI
      RECCO(6) =  COM*SCHI
C---- TRANSFORM TO RECIPROCAL BASE SYSTEM (NB: UB ALWAYS MULTIP.
C     WITH LAMBDA (UBINV DIVIDED BY LAMBDA)
      JUNK = 0
      DO 3048 J=1,3
      DIN(J) = 0.0
      DOUT(J) = 0.0
      DO 3049 K=1,3
      JUNK = JUNK + 1
      DIN(J) = DIN(J) + UBINV(JUNK)*RECCO(K)*WVLNTH
      DOUT(J) = DOUT(J) + UBINV(JUNK)*RECCO(3+K)*WVLNTH 
 3049 CONTINUE
      DOUT(J) = STTH*DOUT(J) - CTTH*DIN(J)
 3048 CONTINUE
C     WRITE(15,9999)IEQPMT,IH,TH,OM,CHI,PHI,DIN,DOUT
      GOTO 20
C
C---- SYNTEX P3 DIFFRACTOMETER
C
C---- FIRST CALCULATE BISECTING ANGLES
C
  400 DO 410 I=1,3                                                             
      RECCO(I)=0                                                                
        DO 410 J=1,3                                                           
        RECCO(I)=RECCO(I)+UB(3*(I-1)+J)*IH(J)                                   
  410 CONTINUE                                                                
C                                                                               
      T1=RECCO(1)*RECCO(1)+RECCO(2)*RECCO(2)                                    
      T2=T1+RECCO(3)*RECCO(3)                                                   
C---- TH IS THETA IN RADIAN.                                                    
      TH=ASIN(0.5*SQRT(T2))                                                     
C---- PHI,CHI IS CLEAR, IS IT NOT                                          
      IF(RECCO(1))4040,4030,4040                                                
 4030 IF(RECCO(2))4040,4035,4040                                                
 4035 PHI=1.5707963                                                            
      GO TO 4045                                                                
 4040 PHI=ATAN2(RECCO(2),RECCO(1))-1.5707963
 4045 CHI=ATAN2(RECCO(3),SQRT(T1))
C     WRITE(15,9998)IH,TH,OM,CHI,PHI,PSI
C9998 FORMAT(1X,3I4,1X,5F8.4)                                             
C---- NOW CALCULATE ANGLES FOR PSI.NE.O
C---- FIRST STORE BISECTING VALUES
      PHIZ = PHI
      CHIZ = CHI
      PSI = PSI*0.017453      
      IF(PSI.EQ.0) GOTO 4046
C---- CALCULATE OMEGA FOR PSI-ROTATION
      IF(ABS(ABS(CHIZ)-1.5707963).LT.1.E-8) GOTO 4046
      OM = ATAN(SIN(PSI)/TAN(CHIZ))
      COSOM = COS(OM)
      GOTO 4047
 4046 COSOM = 1.0
      OM = 0.0
 4047 SINOM = SIN(OM)
C---- CALCULATE TRUE PHI AND CHI
      CHI = ATAN2(SIN(CHIZ)/COSOM,COS(CHIZ)*COS(PSI))
      DELPHI = ATAN2(-SINOM/COS(CHIZ),COSOM*COS(CHI)/COS(CHIZ))
      PHI = PHIZ - DELPHI
C---- PREPARE FOR VECTOR CALCULATION
      SCHI = SIN(CHI)
      CCHI = COS(CHI)
      SPHI = SIN(PHI)
      CPHI = COS(PHI)
      SOM = SIN(TH-OM)
      COM = COS(TH-OM)
      STTH = SIN(2.*TH)
      CTTH = COS(2.*TH)
C---- CALCULATE VECTOR OF BEAM IN AND BEAM OUT IN PHI-SYSTEM
      RECCO(1) =  SOM*CCHI*SPHI - COM*CPHI
      RECCO(2) = -SOM*CCHI*CPHI - COM*SPHI
      RECCO(3) = -SOM*SCHI
      RECCO(4) =  COM*CCHI*SPHI + SOM*CPHI
      RECCO(5) = -COM*CCHI*CPHI + SOM*SPHI
      RECCO(6) = -COM*SCHI
C---- TRANSFORM TO RECIPROCAL BASE SYSTEM (NB: UB ALWAYS MULTIP.
C     WITH LAMBDA (UBINV DIVIDED BY LAMBDA)
      JUNK = 0
      DO 4048 J=1,3
      DIN(J) = 0.0
      DOUT(J) = 0.0
      DO 4049 K=1,3
      JUNK = JUNK + 1
      DIN(J) = DIN(J) + UBINV(JUNK)*RECCO(K)*WVLNTH
      DOUT(J) = DOUT(J) + UBINV(JUNK)*RECCO(3+K)*WVLNTH 
 4049 CONTINUE
      DOUT(J) = STTH*DOUT(J) - CTTH*DIN(J)
 4048 CONTINUE
C     WRITE(15,9999)IEQPMT,IH,TH,OM,CHI,PHI,DIN,DOUT
C                                                                               
C---- NOW TRANSFORMATION FROM RECIPROCAL TO REAL CELL.                          
C     IDENTICAL FOR ALL DIFFRACTOMETER SYSTEMS                       
   20 CALL REZREA(DIN)                                                          
      CALL REZREA(DOUT)                                                         
C                                                                               
C---- CALCULATION OF UNIT VECTOR NORMAL TO THE PLANE DEFINED BY                 
C     VECTORS OF INCIDENT AN REFLECTED BEAM. RESULT IS IN DIN (!).                 
C                                                                               
      CALL ORTHO(DOUT,DIN)                                                      
C                                                                               
C     WRITE(15,9999) IEQPMT,IH,TH,OM,CHI,PHI,DIN,DOUT                                 
C9999 FORMAT(1X,I2,4X,3I4,1X,4F8.4,2X,3F8.4,2X,3F8.4)                                 
      RETURN                                                                    
      END                                                                       
C                                                                               
                       SUBROUTINE REZREA(DIN)                                   
C                                                                               
C     TRANSFORMATION OF THE VECTOR DIN(3)    FROM RECIPROCAL SYSTEM TO          
C                       REAL CRYSTAL SYSTEM.                                    
C                                                                               
C     WRITTEN BY UDO ZUCKER, PROGRAMME TAKES PATTERN FROM XRAY  SYSTEM          
C                                                                               
C                     INPUT ARRAY IS OUTPUT ARRAY                               
C                                                                               
      DIMENSION DIN(3),SD(3)                                                    
C                                                                               
      COMMON/E/DMAX,DMIN,R(6),EXANG(6),EXTRAN(12)                               
C                                                                               
      SD(1)=DIN(1)*R(1)**2+DIN(2)*R(1)*R(2)*R(6)+DIN(3)*R(1)*R(3)*R(5)          
      SD(2)=DIN(1)*R(1)*R(2)*R(6)+DIN(2)*R(2)**2+DIN(3)*R(2)*R(3)*R(4)          
      SD(3)=DIN(1)*R(1)*R(3)*R(5)+DIN(2)*R(2)*R(3)*R(4)+DIN(3)*R(3)**2          
C                                                                               
      DO 20 I=1,3                                                               
   20 DIN(I)=SD(I)                                                              
C                                                                               
      RETURN                                                                    
      END                                                                       
                      SUBROUTINE ORTHO(DIN,DOUT)                                
C                                                                               
C                                                                               
C     COMPUTATION OF UNIT VECTOR NORMAL TO THE PLANE GIVEN BY INCIDENT          
C           AND REFLECTED BEAM OF EXPERIMENT (REAL CRYSTAL SYSTEM).             
C                                                                               
C     BY  UDO ZUCKER, PROGRAMME TAKES PATTERN FROM THE XRAY 80 SYSTEM.          
C                                                                               
C     INPUT VARIABLES                                                           
C     ---------------                                                           
C                                                                               
C     DIN   = ARRAY OF DIM. 3.           VECTOR PARALLEL INCIDENT BEAM.        
C     DOUT  = ARRAY OF DIM. 3.          VECTOR PARALLEL REFLECTED BEAM.
C                                                                               
C     FURTHER PARAMETERS                                                        
C     ------------------                                                        
C                                                                               
C     DR(1)=A-LATTICE CONSTANT   DR(2)=B-LATT.-CON.  DR(3)=C-LATT.CON.          
C     DR(4)=SIN(ALPHA)       DR(5)=SIN(BETA)          DR(6)=SIN(GAMMA)          
C                                                                               
      DIMENSION DIN(3),DOUT(3),DR(6),XN(3),YN(3)                                
C                                                                               
      COMMON/E/DMAX,DMIN,RECELL(6),EXANG(6),EXTRAN(12)                          
      COMMON/FOUR/WINKEL(3),IUB,AG,BG,CG                                        
C                                                                               
      DATA PI/3.1415927/                                                        
C                                                                               
      DR(1)=AG                                                                  
      DR(2)=BG                                                                  
      DR(3)=CG                                                                  
C                                                                               
      DO 100 II=1,3                                                             
  100 DR(3+II) = SIN(WINKEL(II)*PI/180.)                                        
C                                                                               
C---- CALCULATE THE VECTOR PRODUCT  DIN X DOUT IN REAL CRYST. SYSTEM.           
      II=1                                                                      
      JJ=2                                                                      
 5645 K=6-II-JJ                                                                 
      YN(K)= (DIN(II)*DOUT(JJ)-DIN(JJ)*DOUT(II))*DR(II)*DR(JJ)*                 
     &        DR(3+K)/RECELL(K)                                                 
      IF(II.EQ.1) GOTO 5650                                                     
      IF(II.EQ.3) GOTO 5655                                                     
      IF(II.EQ.2) GOTO 5660                                                     
 5650 II=3                                                                      
      JJ=1                                                                      
      GOTO 5645                                                                 
 5655 II=2                                                                      
      JJ=3                                                                      
      GOTO 5645                                                                 
 5660 CONTINUE                                                                  
      DO 106 II=1,3                                                             
  106 XN(II)=YN(II)                                                             
      CALL REZREA(XN)                                                           
C                                                                               
C---- XN IS NOW THE UNNORMALIZED UNIT VECTOR.                                   
C                                                                               
      FNORM=0.0                                                                 
      DO 5665 II=1,3                                                            
 5665 FNORM=FNORM+XN(II)*YN(II)                                                 
      FNORM=SQRT(FNORM)                                                         
C                                                                               
C---- DIN AND DOUT ARE NOW NORMALIZED TO UNIT CRYSTAL AXIS AND DOUT             
C     IN ADDITION TO ONE.                                                       
C                                                                               
      DO 5670 II=1,3                                                            
      DOUT(II)=XN(II)/FNORM*DR(II)                                              
 5670 DIN(II)=DIN(II)*DR(II)                                                    
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
             SUBROUTINE EXTICA(NONPOS,INEUT,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     INEUT  = 0/1  XRAY/NEUTRON RADIATION                                      
C     HELPE  =      COEFFICIENT OF POLARISATION CORRECTION,  CALCULATED         
C                   IN ROUTINE REFINE.                                          
C     HELPF  =      SAME AS HELPE                                               
C                                                                               
C     IMPORTANT VARIABLES                                                       
C     -------------------                                                       
C                                                                               
C     JUNK1 = TYPE OF EXTINCTION                                                
C     IVERT = DISTRIBUTION FUNCTION SELECTOR                                    
C     FCALC = IS HERE FCALC*FCALC                                               
C                                                                               
C                                                                               
      DIMENSION VECTOR(3),TERM(3,3),QTQ(35),IQTQ(35),IBDF(2)                    
C                                                                               
      COMMON/D/ACALC,AIND(60),ANOM(2,20),BCALC,BIND(60),DADB(6,60),DADO(        
     & 60),DADX(3,60),DBDB(6,60),DBDO(60),DBDX(3,60),P(3),IH(3),ISCAT(2,        
     & 60),ISCL,MFSEL,OCCA(60),OCCUP(60), SCALE(10), SCAT(20), SITE(60),        
     & DADC(10,60),DBDC(10,60),DADD(15,60), DBDD(15,60), AMINUS, BMINUS,        
     & DRDE(7),NDEGN,VORFAK,IPAREX(7),DAALFA(60),DBALFA(60)                     
      COMMON/G/BUFF(1023),FCURVE(9,20),IFMT(18),IWT,NSCAT,BETFAC,WVLNTH         
     &,WTFAC                                                                    
      COMMON/E/DMAX,DMIN,RECELL(6),EXANG(6),EXTRAN(12)                          
      COMMON/F/DELTA,FCALC,FOBS,FOMAX,FOMIN,IREJ,MREF,MREJ,NREJ,NUMBER(3        
     & 6), RDEN(36),RHO,RHOMN,RHOMX,RNUM(36),WT, EXBETA,EXX,EXY, GAMMAP,        
     & FISGN,IFSGN,EXTIN(6),KEXTRF(7),DELMAX,IANIGC,GAMMAS,IZACH,RNEUEX(        
     & 7),NEUEXT,NDISTR,NIEXT,WTA(3),WTB(3),RHODIV                              
      COMMON/COPPEN/YEXT,BDF(16),JUNK6                                          
C                                                                               
      EQUIVALENCE (QTQ,IQTQ),(BDF,IBDF)                                         
C                                                                               
C                                                                               
      JUNK1 = IBDF(1)+1                                                         
      IVERT = IBDF(2)                                                           
C                                                                               
C---- THE VALUE GAMMAS IS IN BINARY STORAGE CALLED EXBETA.  YEXT IS EX -        
C     TINCTION FACTOR.                                                          
C                                                                               
      GAMMS=EXBETA                                                              
      YEXT=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)*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                                                                               
      GAMMS=EXBETA                                                              
C                                                                               
C---- POLARISATIONCORRECTION                                                    
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(5779,5719,5716,5719,5690,5690,5690),JUNK1                            
C                                                                               
C---- ANISOTROPIC, CALCULATE VECTOR(TRANSPOSE)*TENSOR*VECTOR.                   
C                                                                               
 5690 DO 5695 I=1,6                                                             
 5695 QTQ(I)=BDF(4+I)                                                           
      IF(JUNK1.NE.6) GOTO 5703                                                  
      DO 5700 I=1,6                                                             
 5700 QTQ(I)=BDF(10+I)*WVLNTH**2                                                
 5703 DO 5706 II=1,3                                                            
      VECTOR(II)=EXANG(3+II)                                                    
      IF(JUNK1.EQ.6)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.GT.0.) GOTO 5710                                                  
C                                                                               
C---- RESULT IS INPHYSICALLY.                                                   
      NONPOS=2                                                                  
      H=0.                                                                      
      GOTO 5713                                                                 
C                                                                               
C---- RESULT IS MEANINGFUL                                                      
C                                                                               
 5710 H=1.0/SQRT(PROD)                                                          
 5713 CONTINUE                                                                  
C                                                                               
C---- HERE JUNK IS GREATER  4.                                                  
C                                                                               
      GOTO(5901,5719,5716,5719,5719,5716,5719),JUNK1                            
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 AS=AG                                                                     
      BS=BG                                                                     
      CS=2.12                                                                   
      IF(IVERT.EQ.0) GOTO 5723                                                  
      AS=AL                                                                     
      BS=BL                                                                     
      CS=2.0                                                                    
 5723 CONTINUE                                                                  
C                                                                               
C---- CALCULATE DERIVATIVES                                                     
C                                                                               
      GOTO(5901,5729,5726,5736,5729,5726,5736),JUNK1                            
C                                                                               
C---- TYPE 2  DERIVATIVE  DBRHO IN QTQ(10)                                      
C                                                                               
 5726 RO=BDF(4)/WVLNTH                                                          
      IF(JUNK1.EQ.6) 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(JUNK1.EQ.5) 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)/WVLNTH                                                          
      GPAR=BDF(3)                                                               
      IF(JUNK1.EQ.4) 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(JUNK1.EQ.7.OR.JUNK1.EQ.4)     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 TRETED.                                    
C                                                                               
      QTQ(12)=0.0                                                               
      IF(JUNK1.EQ.4.OR.JUNK1.EQ.7) QTQ(12)=ABLEIT(AP,BP,CP,XP)                  
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,5756,5753                                                   
C                                                                               
C---- NEUTRON                                                                   
C                                                                               
 5753 YEXT=YP*YS                                                                
      QTQ(16)=YEXT*(1.+XP*QTQ(12)+XS*QTQ(13)*(1.+XP*QTQ(12)))                   
      QTQ(17)=YEXT*YP*QTQ(13)                                                   
      QTQ(18)=YEXT*(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 )                            
      QTQ(14)=0.0                                                               
      IF(JUNK1.EQ.4.OR.JUNK1.EQ.7) QTQ(14)=ABLEIT(AP,BP,CP,XPP)                 
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                                                                               
      YEXT=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                                                                               
C---- REFINEMENT ON F**2                                                        
C                                                                               
      QTQ(19)=FCALC**2                                                          
C                                                                               
C---- REFINEMENT ON F                                                           
C                                                                               
      IF(MREF.EQ.0) QTQ(19)=0.5*((SQRT(FCALC))**3)/(SQRT(YEXT))                 
C                                                                               
C---- NOW COMPUTATION IS SAME FOR BOTH CASES.                                   
C                                                                               
      GOTO(5901,5763,5763,5763,5766,5766,5766),JUNK1                            
C                                                                               
C---- ISOTROPIC CASE ; DERIVATIVE IN QTQ(20) (+IN QTQ(21) IN G.CASE)            
C                                                                               
 5763 IF(JUNK1.NE.2) QTQ(20)=QTQ(19)*QTQ(18)/WVLNTH                             
      IF(JUNK1.EQ.2) QTQ(20)=QTQ(19)*QTQ(17)*GAMMS*QTQ(11)                      
      IF(JUNK1.EQ.4) QTQ(21)=QTQ(19)*QTQ(17)*GAMMS*QTQ(11)                      
      GOTO 5773                                                                 
C                                                                               
C---- ANISOTROPIC CASE. DERIVATIVESIN  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(JUNK1.NE.6)QTQ(19+K)=QTQ(19)*GAMMS*QTQ(17)*QTQ(11)*TERM(I,J)*          
     &QTQ(27)                                                                   
      IF(JUNK1.EQ.6)QTQ(19+K)=QTQ(19)*QTQ(18)*TERM(I,J)*QTQ(27)*WVLNTH**        
     &2                                                                         
CCC   QTQ(20)=QTQ(20)+QTQ(21)-0.5*QTQ(23)                                       
 5769 CONTINUE                                                                  
      IF(JUNK1.EQ.7) QTQ(26)=QTQ(19)*QTQ(18)/WVLNTH                             
 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.GT.0) GOTO 5776                                                   
C---- REFINEMENT ON F   >>>> VORFAKTOR=  OMEGA1 / FCALC*SQRT(YEXT)              
C                                                                               
      VORFAK= QTQ(16)/SQRT(FCALC*YEXT)                                          
      GOTO 5779                                                                 
C                                                                               
C *** REFINEMENT ON F**2 >>>> VORFAKTOR= OMEGA1*2                               
C                                                                               
 5776 VORFAK=QTQ(16)*2                                                          
 5779 CONTINUE                                                                  
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 (5999,5956,5956,5956,5963,5963,5963),JUNK1                           
C                                                                               
C---- ISOTROPIC CASES                                                           
C                                                                               
 5956 JUNK3 = 1                                                                 
      IF(JUNK1.EQ.4) JUNK3=2                                                    
      DO 5959 II=1,JUNK3                                                        
 5959 DRDE(II)=QTQ(19+II)                                                       
      GOTO 5999                                                                 
C                                                                               
C---- ANISOTROPIC CASES                                                         
C                                                                               
 5963 JUNK3=6                                                                   
      IF(JUNK1.EQ.7) JUNK3=7                                                    
      DO 5966 II=1,JUNK3                                                        
 5966 DRDE(II)=QTQ(19+II)                                                       
 5999 RETURN                                                                    
      END                                                                       
C                                                                               
      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                                                                       
C                                                                               
      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                                                                       
C                                                                               
                            SUBROUTINE WEIGHT                                   
C                                                                               
C          WRITTEN BY UDO ZUCKER  MPI-FKF STUTTGART  OCTOBER 1979               
C                                                                               
C     THIS ROUTINE HAS TWO ENTRY-POINTS. THE FIRST ONE IS FOR  COMPUTA-         
C     TION OF WEIGHTS.      THE   SECOND IS FOR CALCULATION WHETHER   A         
C     REFLECTION  IS INCLUDED IN MATRIX-INVERSION OR NOT. SEVERAL CON -         
C     DITIONS ARE AVAILABLE.                                                    
C                                                                               
      COMMON/LIMIT/UNTEN,OBEN,IREJEC                                            
      COMMON/F/DELTA,FCALC,FOBS,FOMAX,FOMIN,IREJ,MREF,MREJ,NREJ,NUMBER(3        
     & 6),RDEN(36),RHO,RHOMN,RHOMX,RNUM(36),WT,EXBETA,EXX,EXY,GAMMAP            
     &,FISGN,IFSGN,EXTIN(6),KEXTRF(7),DELMAX,IANIGC,GAMMAS,IZACH,RNEUEX(        
     &7),NEUEXT,NDISTR,NIEXT,WTA(3),WTB(3),RHODIV                               
C                                                                               
C---- SPECIAL WEIGHT ROUTINE ADDED BY W.F.KUHS MAY 82                           
C                                                                               
C                                                                               
C---- FIRST APPLY REES FORMALISM TO SIGMA(F**2)                                 
C                                                                               
      IF(FOBS*FOBS-WT)20,20,10                                                  
   10 WT = 0.5*WT/FOBS                                                          
      GOTO 30                                                                   
   20 WT = SQRT(WT)*0.5                                                         
C                                                                               
C---- APPLY DIFFERENT WEIGTHS DEPENDING ON SINTH/LAM                            
C                                                                               
   30 ROOT = SQRT(RHO)                                                          
      IF(ROOT-RHODIV) 50,50,40                                                  
   40 WT = WT*(WTB(1)+WTB(2)*ROOT+WTB(3)*RHO)                                   
      GOTO 60                                                                   
   50 WT = WT*(WTA(1)+WTA(2)*ROOT+WTA(3)*RHO)                                   
   60 RETURN                                                                    
C                                                                               
C                                                                               
C **********************************************************************        
      ENTRY REJECT                                                              
C **********************************************************************        
C                                                                               
C     RHO=(SIN(THETA)/LAMBDA)**2                                                
      ROOT=SQRT(RHO)                                                            
C                                                                               
      GOTO(1,2,3,4),IREJEC                                                      
C                                                                               
C---- INCLUDE ALL REFLECTION WHICH ARE NOT  SIGNED AS *L-THAN* AND WHICH        
C     ARE LAYING IN THE SIN(THETA)/LAMDA INTERVAL  FROM  UNTEN  TO OBEN.        
C                                                                               
    1 IF(UNTEN-ROOT.GT.0..OR.ROOT-OBEN.GT.0..OR.MREJ.EQ.1)  NREJ=1              
      GOTO 5                                                                    
C                                                                               
C---- INCLUDE THOSE *LTHAN-REFLECTIONS*,WHOSE FCALC IS GREATER FOBS.            
C                                                                               
    2 IF(MREJ.EQ.0)GOTO 5                                                       
      IF((FCALC-FOBS).LT.0.) NREJ=1                                             
      GOTO 5                                                                    
C                                                                               
C---- INCLUDE ALL REFLECTIONS, WHICH ARE LAYING IN THE SIN(THETA)/LAMBDA        
C     INTERVAL  FROM  UNTEN  TO  OBEN   AND FOR WHICH FCALC GT FOBS  IS         
C     VALID.                                                                    
C                                                                               
    3 IF(UNTEN-ROOT.GT.0..OR.ROOT-OBEN.GT.0.)  NREJ=1                           
      IF(NREJ.EQ.1.OR.MREJ.EQ.0) GOTO 5                                         
      IF((FCALC-FOBS).LT.0.) NREJ=1                                             
      GOTO 5                                                                    
C                                                                               
C---- DUMMY-ROUTINE FOR SPECIAL REQUIREMENTS                                    
    4 GOTO 5                                                                    
C                                                                               
    5 RETURN                                                                    
      END                                                                       
C                                                                               
      SUBROUTINE GERMIT (HPR,QTRANS,ALF,HHJ,HHHH,ISOT,KSYM,T1,T2,TY)            
C                                                                               
C     ROUTINE FOR CALCULATION  OF POLYNOMS G(JKL) AND G(JKLM), WHICH            
C     ARE RELATED TO THE HERMITE  POLYNOMS.   IF ALPHA = 0 ONE  GETS            
C     HERMITE POLYNOMS. FOR JOKE  THIS POLYNOMS ARE CALLED *GERMITS*            
C                                                                               
C     H1,H2,H3 ARE MILLER INDICES,  MULTIPLES OF PI ARE INCLUDED  IN            
C            THE COEFFICIENTS FAK1, FAK2 IN ROUTINE SFAC.                       
C                                                                               
C                       WRITTEN BY UDO ZUCKER                                   
C                                                                               
C     INPUT VARIABLES                                                           
C     ---------------                                                           
C                                                                               
C     HPR   =  INPUT  ARRAY CONTAINING H,K, L INDICES OF TRANSFORMED            
C              REFLECTION.                                                      
C     QTRANS=  INPUT  ARRAY CONTAINING  BETA  MATRIX OF  TRANSFORMED            
C              REFLECTION.                                                      
C     ALF   =  ALPHA PARAMETER OF CURRENT ATOM.                                 
C     HHJ   =  OUTPUT ARRAY CONTAIN. HERMITE POLYN.  OF THIRD ORDER.            
C     HHHH  =  OUTPUT ARRAY CONTAIN. HERMITE POLYN. OF FOURTH ORDER.            
C     ISOT  =  TEMPERATURE FACTOR TYPE  INDICATOR OF  CURRENT  ATOM.            
C     KSYM  =  NUMBER OF CURRENT SYMMETRY OPERATION.                            
C     T1,T2                                                                     
C     AND TY=  TEMPORARY FIELDS CONTAINING INFORMATION ABOUT  DERI-             
C              VATIVES OF HERMITE POLYNOMS WTH TO ALPHA  ( SEE SUB-             
C              ROUTINE ABALFA).                                                 
C                                                                               
      DIMENSION HPR(3),Q(3,3),QTRANS(6,24),HHJ(10),HHHH(15), T1(15),            
     &T2(15),TY(10)                                                             
C                                                                               
      H1  = HPR(1)                                                              
      H2  = HPR(2)                                                              
      H3  = HPR(3)                                                              
      H11 = H1*H1                                                               
      H22 = H2*H2                                                               
      H33 = H3*H3                                                               
C                                                                               
C---- COPY QTRANS IN Q FOR EASIER UNDERSTANDING.                                
C                                                                               
      DO 100 K=1,3                                                              
  100 Q(K,K)=QTRANS(K,KSYM)                                                     
      Q(1,2)=QTRANS(4,KSYM)                                                     
      Q(1,3)=QTRANS(5,KSYM)                                                     
      Q(2,3)=QTRANS(6,KSYM)                                                     
C                                                                               
      ALFQ=ALF*ALF                                                              
C                                                                               
C---- GERMIT-POLYNOMS OF THIRD ORDER.                                           
C                                                                               
      HHJ(1) = H11*H1*1.E-3                                                     
      HHJ(2) = H22*H2*1.E-3                                                     
      HHJ(3) = H33*H3*1.E-3                                                     
      HHJ(4) = H11*H2*3.E-3                                                     
      HHJ(5) = H22*H1*3.E-3                                                     
      HHJ(6) = H11*H3*3.E-3                                                     
      HHJ(7) = H1*H33*3.E-3                                                     
      HHJ(8) = H22*H3*3.E-3                                                     
      HHJ(9) = H2*H33*3.E-3                                                     
      HHJ(10)= H1*H2*H3*6.E-3                                                   
C                                                                               
      TY (1) =  -3.E-3 * H1*Q(1,1)                                              
      TY (2) =  -3.E-3 * H2*Q(2,2)                                              
      TY (3) =  -3.E-3 * H3*Q(3,3)                                              
      TY (4) = (-2.*H1*Q(1,2)-H2*Q(1,1))*3.E-3                                  
      TY (5) = (-H1*Q(2,2)-2.*H2*Q(1,2))*3.E-3                                  
      TY (6) = (-2.*H1*Q(1,3)-H3*Q(1,1))*3.E-3                                  
      TY (7) = (-H1*Q(3,3)-2.*H3*Q(1,3))*3.E-3                                  
      TY (8) = (-2.*H2*Q(2,3)-H3*Q(2,2))*3.E-3                                  
      TY (9) = (-H2*Q(3,3)-2.*H3*Q(2,3))*3.E-3                                  
      TY (10)= (-H1*Q(2,3)-H2*Q(1,3)-H3*Q(1,2))*6.E-3                           
C                                                                               
      IF(ABS(ALF).LT.1.E-8) GOTO 3                                              
C                                                                               
      DO 5 I = 1,10                                                             
    5 HHJ(I) = HHJ(I) + TY(I)*ALF                                               
C                                                                               
    3 IF(ISOT-2)1,1,2                                                           
C                                                                               
C---- *GERMITS* OF FOURTH ORDER.                                                
C                                                                               
    2 HHHH(1)  = 1. E-4 * H11*H11                                               
      HHHH(2)  = 1. E-4 * H22*H22                                               
      HHHH(3)  = 1. E-4 * H33*H33                                               
      HHHH(4)  = 4. E-4 * H11*H1*H2                                             
      HHHH(5)  = 4. E-4 * H11*H1*H3                                             
      HHHH(6)  = 6. E-4 * H11*H22                                               
      HHHH(7)  = 12.E-4 * H11*H2*H3                                             
      HHHH(8)  = 6. E-4 * H11*H33                                               
      HHHH(9)  = 4. E-4 * H22*H2*H1                                             
      HHHH(10) = 12.E-4 * H22*H3*H1                                             
      HHHH(11) = 12.E-4 * H1*H2*H33                                             
      HHHH(12) = 4. E-4 * H33*H3*H1                                             
      HHHH(13) = 4. E-4 * H22*H2*H3                                             
      HHHH(14) = 6. E-4 * H33*H22                                               
      HHHH(15) = 4. E-4 * H33*H3*H2                                             
C                                                                               
                                                                                
      T1  (1) = -6. E-4 * H11*Q(1,1)                                            
      T2  (1) = +3. E-4 * Q(1,1)*Q(1,1)                                         
      T1  (2) = -6. E-4 * H22*Q(2,2)                                            
      T2  (2) =  3. E-4 * Q(2,2)*Q(2,2)                                         
      T1  (3) = -6. E-4 * H33*Q(3,3)                                            
      T2  (3) =  3. E-4 * Q(3,3)*Q(3,3)                                         
      T1  (4) = -12.E-4 * (H11*Q(1,2) + H1*H2*Q(1,1))                           
      T2  (4) =  12.E-4 * Q(1,1)*Q(1,2)                                         
      T1  (5) = -12.E-4 * (H11*Q(1,3) + H1*H3*Q(1,1))                           
      T2  (5) =  12.E-4 * Q(1,1)*Q(1,3)                                         
      T1  (6) = -6. E-4 * (H11*Q(2,2) + 4.*H1*H2*Q(1,2) + H22*Q(1,1))           
      T2  (6) =  6. E-4 * (Q(1,1)*Q(2,2)+2.*Q(1,2)*Q(1,2))                      
      T1  (7) = -12.E-4 * (H11*Q(2,3) + 2.*H1*H2*Q(1,3) + 2.*H1*H3*             
     &                     Q(1,2) + H2*H3*Q(1,1))                               
      T2  (7) =  12.E-4 * (Q(1,1)*Q(2,3)+2.*Q(1,2)*Q(1,3))                      
      T1  (8) = -6. E-4 * (H11*Q(3,3)+4.*H1*H3*Q(1,3)+H33*Q(1,1))               
      T2  (8) =  6. E-4 * (Q(1,1)*Q(3,3)+2.*Q(1,3)*Q(1,3))                      
      T1  (9) = -12.E-4 * (H22*Q(1,2) + H2*H1*Q(2,2))                           
      T2  (9) =  12.E-4 * Q(2,2)*Q(1,2)                                         
      T1  (10)= -12.E-4 * (H22*Q(1,3)+2.*H2*H3*Q(1,2)+2.*H2*H1*Q(2,3)           
     &                    +H3*H1*Q(2,2))                                        
      T2  (10)=  12.E-4 * (Q(2,2)*Q(1,3) + 2.*Q(2,3)*Q(1,2))                    
      T1  (11)= -12.E-4 * (H1*H2*Q(3,3)+H2*H3*Q(1,3)*2.+H33*Q(1,2)              
     &                     +2.*H1*H3*Q(2,3))                                    
      T2  (11)=  12.E-4 * (Q(1,2)*Q(3,3)+2.*Q(1,3)*Q(2,3))                      
      T1  (12)= -12.E-4 * (H33*Q(1,3) + H1*H3*Q(3,3))                           
      T2  (12)=  12.E-4 * Q(3,3)*Q(1,3)                                         
      T1  (13)= -12.E-4 * (H22*Q(2,3) + H2*H3*Q(2,2))                           
      T2  (13)=  12.E-4 *  Q(2,2)*Q(2,3)                                        
      T1  (14)= -6. E-4 * (H33*Q(2,2)+4.*H3*H2*Q(2,3)+H22*Q(3,3))               
      T2  (14)=  6. E-4 * (Q(3,3)*Q(2,2)+2.*Q(2,3)*Q(2,3))                      
      T1  (15)= -12.E-4 * (H33*Q(2,3) + H2*H3*Q(3,3))                           
      T2  (15)=  12.E-4 * Q(3,3)*Q(2,3)                                         
C                                                                               
      IF(ABS(ALF).LT.1.E-8) GOTO 1                                              
C                                                                               
      DO 10 I = 1,15                                                            
   10 HHHH(I) = HHHH(I) + ALF*T1(I) + ALFQ*T2(I)                                
C                                                                               
    1 RETURN                                                                    
      END                                                                       
C                                                                               
                      SUBROUTINE UMKEHR(BETA,I,Q)                               
C                                                                               
C     ROUTINE  INVERTS  THE ANISOTROPIC TEEMPERATURE FACTOR BETA OF             
C     THE I TH ATOM AND  STORES  THE RESULT IN USUAL MANNER INTO Q.             
C                                                                               
C                          BY UDO ZUCKER                                        
C                                                                               
      DIMENSION BETA(6,60),Q(3,3)                                               
C                                                                               
      Q(1,1) = BETA(2,I)*BETA(3,I)-BETA(6,I)**2                                 
      Q(1,2) = -BETA(4,I)*BETA(3,I)+BETA(5,I)*BETA(6,I)                         
      Q(1,3) = BETA(4,I)*BETA(6,I)-BETA(5,I)*BETA(2,I)                          
      QDET   = BETA(1,I)*Q(1,1)+BETA(4,I)*Q(1,2)+BETA(5,I)*Q(1,3)               
      Q(1,1) = Q(1,1)/QDET                                                      
      Q(1,2) = Q(1,2)/QDET                                                      
      Q(1,3) = Q(1,3)/QDET                                                      
      Q(2,2) = (BETA(1,I)*BETA(3,I)-BETA(5,I)**2)/QDET                          
      Q(2,3) = (-BETA(1,I)*BETA(6,I)+BETA(5,I)*BETA(4,I))/QDET                  
      Q(3,3) = (BETA(1,I)*BETA(2,I)-BETA(4,I)**2)/QDET                          
      Q(2,1) = Q(1,2)                                                           
      Q(3,1)=Q(1,3)                                                             
      Q(3,2)=Q(2,3)                                                             
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
                SUBROUTINE ABALFA(DHH4,ALF,T1,T2,ISOT)                          
C                                                                               
C     COMPUTATION OF THE DERIVATIVES OF THE *GERMIT* POLYNOMS WORTH             
C                          TO ALPHA.                                            
C                                                                               
C           WRITTEN BY  UDO  ZUCKER (MPI/FKF) JULY 79                           
C           REVISED  AND  RESEQUENCED  DECEMBER  1980                           
C                                                                               
      DIMENSION DHH4(15),T1(15),T2(15)                                          
C                                                                               
      TWOALF = ALF*2.                                                           
C                                                                               
C---- DERIVATIVES OF *GERMITS* OF FOURTH ORDER WTH TO ALPHA.                    
C                                                                               
      IF(ISOT-2) 1,1,2                                                          
C                                                                               
    2 DO 3 I = 1,15                                                             
    3 DHH4(I)= T1(I) + TWOALF*T2(I)                                             
C                                                                               
    1 RETURN                                                                    
      END                                                                       
                                                                                
