C     ==================================================================        
C     !!                                                              !!        
C     !!    DDDDDDD   IIII   CCCCCC   VV    VV   OOOOOO   LL          !!        
C     !!    DD    DD   II   CC    CC  VV    VV  OO    OO  LL          !!        
C     !!    DD    DD   II   CC        VV    VV  OO    OO  LL          !!        
C     !!    DD    DD   II   CC        VV    VV  OO    OO  LL          !!        
C     !!    DD    DD   II   CC         VV  VV   OO    OO  LL          !!        
C     !!    DD    DD   II   CC    CC    VVVV    OO    OO  LL    LL    !!        
C     !!    DDDDDDD   IIII   CCCCCC      VV      OOOOOO   LLLLLLLL    !!        
C     !!                                                              !!        
C     ==================================================================        
C                                                                               
C     TRIAL AND ERROR METHOD FOR THE AUTOMATIC INDEXING OF POWDER               
C     PATTERNS - VARIATION OF PARAMETERS BY SUCCESSIVE DICHOTOMY.               
C     THIS VERSION INCLUDES PARTITION OF VOLUME SPACE                           
C     VERSION : 30-09-1982                                                      
C     THE FOLLOWING CRYSTALLINE SYSTEMS CAN BE EXAMINED :                       
C     CUBIC, QUADRATIC, HEXAGONAL, ORTHORHOMBIC AND MONOCLINIC SYSTEMS.         
C                                                                               
C     FOR FURTHER INFORMATION CONTACT                                           
C               D. LOUER                                                        
C               LABORATOIRE DE CRISTALLOCHIMIE                                  
C               UNIVERSITE DE RENNES I                                          
C               CAMPUS DE BEAULIEU                                              
C               35042 RENNES CEDEX    -FRANCE-                                  
C               TEL.  16-99-364815-(POSTE 2258)                                 
C                                                                               
C                                                                               
C      ---------------------------------------------------------------          
C      !            G E N E R A L   I N S T R U C T I O N S          !          
C      ---------------------------------------------------------------          
C                                                                               
C      DICVOL PROPOSES SOLUTIONS; THE USER DISPOSES OF THEM.                    
C      THE AUTHORS DISCLAIM RESPONSABILITY FOR ALL USE OF THIS                  
C      PROGRAM, WHETHER GOOD, BAD OR INDIFFERENT.                               
C                                                                               
C      DATA CARDS                                                               
C      ----------                                                               
C                                                                               
C  CARD 1  TITLE                                       FORMAT(16A4)             
C                                                                               
C                                                                               
C                                                                               
C  CARD 2  N,ITYPE,JC,JT,JH,JO,JM                      FORMAT(7I2)              
C                                                                               
C          N     (COL. 1- 2)  NUMBER OF LINES (E.G. 20).                        
C          ITYPE (COL.    4)  SPACING DATA TYPE.                                
C                         =1  THETA BRAGG IN DEGREES.                           
C                         =2  2-THETA ANGLE IN DEGREES.                         
C                         =3  D-SPACING IN ANGSTROMS.                           
C                         =4  Q SPECIFIED IN Q-UNITS AS E+04/D**2.              
C          JC    (COL.    6)                                                    
C                         =0  CUBIC SYSTEM IS NOT TESTED.                       
C                         =1  CUBIC SYSTEM IS TESTED.                           
C          JT    (COL.    8)                                                    
C                         =0  TETRAGONAL SYSTEM IS NOT TESTED.                  
C                         =1  TETRAGONAL SYSTEM IS TESTED.                      
C          JH    (COL.   10)                                                    
C                         =0  HEXAGONAL SYSTEM IS NOT TESTED.                   
C                         =1  HEXAGONAL SYSTEM IS TESTED.                       
C          JO    (COL.   12)                                                    
C                         =0  ORTHORHOMBIC SYSTEM IS NOT TESTED.                
C                         =1  ORTHORHOMBIC SYSTEM IS TESTED.                    
C          JM    (COL.   14)                                                    
C                         =0  MONOCLINIC SYSTEM IS NOT TESTED.                  
C                         =1  MONOCLINIC SYSTEM IS TESTED.                      
C                                                                               
C                                                                               
C  CARD 3  AMAX,BMAX,CMAX,VOLMIN,VOLMAX,BEMIN,BEMAX    FORMAT(7F10.2)           
C                                                                               
C          AMAX  (COL. 1-10)  MAXIMUM VALUE OF UNIT CELL DIMENSION A            
C                                        IN ANGSTROMS.                          
C          BMAX  (COL.11-20)  MAXIMUM VALUE OF UNIT CELL DIMENSION B            
C                                        IN ANGSTROMS.                          
C          CMAX  (COL.21-30)  MAXIMUM VALUE OF UNIT CELL DIMENSION C            
C                                        IN ANGSTROMS.                          
C          VOLMIN(COL.31-40)  MINIMUM VOLUME FOR TRIAL UNIT CELLS               
C                                        IN ANGSTROMS**3.                       
C          VOLMAX(COL.41-50)  MAXIMUM VOLUME FOR TRIAL UNIT CELLS               
C                                        IN ANGSTROMS**3.                       
C          BEMIN (COL.36-42)  MINIMUM ANGLE FOR UNIT CELL IN DEGREES            
C                                    (DEFAULT= 90 DEGREES).                     
C          BEMAX (COL.43-49)  MAXIMUM ANGLE FOR UNIT CELL IN DEGREES            
C                                    (DEFAULT=130 DEGREES).                     
C                                                                               
C                                                                               
C  CARD 4  WAVE,POIMOL,DENS,DELDEN                     FORMAT(4F10.5)           
C                                                                               
C          WAVE  (COL. 1-10)  WAVELENGTH IN ANGSTROMS.                          
C          POIMOL(COL.11-20)  MOLECULAR WEIGHT OF ONE FORMULA UNIT IN           
C                             A.M.U.  (DEFAULT =0.0 IF FORMULA WEIGHT           
C                             NOT KNOWN).                                       
C          DENS  (COL.21-30)  MEASURED DENSITY IN G.CM(-3)                      
C                             (DEFAULT =0.0 IF DENSITY NOT KNOWN).              
C          DELDEN(COL.31-41)  ABSOLUTE ERROR IN MEASURED DENSITY.               
C                                                                               
C                                                                               
C  CARD 5 TO 5+N  D(I),EPSIL(I)                        FORMAT(2F10.4)           
C                                                                               
C          (ONE FOR EACH OBSERVED LINE, UP TO N)                                
C          D(I)  (COL. 1-10)  VALUE DESCRIBING THE OBSERVED POSITION            
C                             OF THIS LINE ACCORDING TO 'ITYPE'.                
C          EPSIL (COL.11-20)  ABSOLUTE ERROR IN 'D(I)', ACCORDING TO            
C                             'ITYPE'.                                          
C                                                                               
C  NOTE:                                                                        
C          IF ITYPE=1,2,4 THE VALUES OF 'D(I)' AND 'EPSIL(I)' MUST              
C          BE LISTED IN INCREASING ORDER.                                       
C          IF ITYPE=3 THEY MUST BE IN DECREASING ORDER.                         
C                                                                               
C  REFERENCES:                                                                  
C     1.-  LOUER, D. AND LOUER, M., METHODE D'ESSAIS ET ERREURS POUR            
C          L'INDEXATION AUTOMATIQUE DES DIAGRAMMES DE POUDRE,                   
C          J. APPL. CRYST. 5,271-275(1972).                                     
C     2.-  LOUER, D., SUCCESSIVE DICHOTOMY METHOD FOR INDEXING POWDER           
C          PATTERNS, SYMPOSIUM ON ACCURACY IN POWDER DIFFRACTION, N.B.S.        
C          GAITHERSBURG, MD., U.S.A., JUNE 11-15,1979.                          
C     3.-  LOUER,D. AND VARGAS, R., INDEXATION AUTOMATIQUE DES                  
C          DIAGRAMMES DE POUDRE PAR DICHOTOMIES SUCCESSIVES,                    
C          J. APPL. CRYST. 15,542-545,(1982).                                   
C     4.-  DE WOLFF, P. M., A SIMPLIFIED CRITERION FOR THE RELIABILITY          
C          OF A POWDER PATTERN INDEXING, J. APPL. CRYST. 5,108-113(1968)        
C     5.-  SMITH, G. S. AND SNYDER, R. L., F(N): A CRITERION FOR                
C          RATING POWDER DIFFRACTION PATTERNS AND EVALUATING THE RELIA-         
C          BILITY OF POWDER-PATTERN INDEXING, J. APPL. CRYST. 12,60-65          
C          (1979).                                                              
C                                                                               
C                                                                               
C           ----------          ------------------          ----------          
C           !        !=========>!     CUBIQUE    !=========>!        !          
C           !        !          ------------------          !        !          
C           !        !                   V                  !        !          
C           !        !              ----------              !        !          
C           !        !              ! CUDHKL !              !        !          
C           !        !              ----------              !    A   !          
C           !        !                                      !    F   !          
C           !        !                                      !    F   !          
C           !        !                                      !    I   !          
C           !        !          ------------------          !    N   !          
C           !        !=========>!   TETRAGONAL   !=========>!    E   !          
C           !    P   !          ------------------          !    M   !          
C           !    R   !                   V                  !    E   !          
C           !    O   !              ----------              !    N   !          
C           !    G   !              ! TEDHKL !              !    T   !          
C           !    R   !              ----------              !        !          
C           !    A   !                                      !        !          
C           !    M   !                                      !    D   !          
C           !    M   !                                      !    E   !          
C           !    E   !          ------------------          !    S   !          
C           !        !=========>!    HEXAGONAL   !=========>!        !          
C           !        !          ------------------          !        !          
C           !    P   !                   V                  !    P   !          
C           !    R   !              ----------              !    A   !          
C           !    I   !              ! TEDHKL !              !    R   !          
C           !    N   !              ----------              !    A   !          
C           !    C   !                                      !    M   !          
C           !    I   !                                      !    E   !          
C           !    P   !                                      !    T   !          
C           !    A   !          ------------------          !    R   !          
C           !    L   !=========>! ORTHORHOMBIQUE !=========>!    E   !          
C           !        !          ------------------          !    S   !          
C           !        !                   V                  !        !          
C           !        !              ----------              !        !          
C           !        !              ! ORDHKL !              !        !          
C           !        !              ----------              !        !          
C           !        !                                      !        !          
C           !        !                                      !        !          
C           !        !                                      !        !          
C           !        !          ------------------          !        !          
C           !        !=========>!  MONOCLINIQUE  !=========>!        !          
C           !        !          ------------------          ----------          
C           !        !                   V                       V              
C           !        !              ----------              ----------          
C           !        !              ! MODHKL !              ! PASAJE !          
C           ----------              ----------              ----------          
C     ------------------------------------------------------------------        
C     !              P R O G R A M M E   P R I N C I P A L             !        
C     ------------------------------------------------------------------        
      DIMENSION IDENT(16)                                                       
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      CHARACTER*40 FILEX
C     CHANGES GGJJ MAY 13,1986
C     MORE CHANGES 12/87 RLS
      WRITE(6,5001)
5001  FORMAT(' VAX VERSION DEC 1987---LOUER CODE OF SEPT 30, 1982'/
     1     10X,'MRL/PSU MODS'/
     1     10X,'BY GERALD G. JOHNSON, JR,'/
     1     10X,'AU MODS BY ROBERT L. SNYDER'/
     2     10X,'TO LOUER"S DICVOL INDEXING PROGRAM'//)
C5003  WRITE(6,5002)
C      FILEX='      '
C5002  FORMAT(' PLEASE ENTER THE NAME OF THE INPUT DATA FILE.'/
C     1    10X,'OR THE WORD "HELP" TO SEE A LIST OF *.LOU FILES'/
C     2    10X,'OR THE WORD "NOFILE" TO ABORT     ',$)
C      READ(5,5004,END=9977,ERR=9977) FILEX
C5004  FORMAT(A)
C      IKEY=STR$UPCASE(FILEX,FILEX)
C      IF(FILEX.EQ.'      ') FILEX='HELP'
C      IF(FILEX.EQ.'NOFILE') STOP ' JOB ABORTED--START AGAIN'
C      IF(FILEX.EQ.'HELP')   THEN
C      IKEY=LIB$SPAWN('$ DIR *.LOU')
C      GO TO 5003
C      ENDIF
C      OPEN(UNIT=5,READONLY,FILE=FILEX,ERR=5003,STATUS='OLD',
C     1      DEFAULTFILE='.LOU',SHARED)
C      OPEN(UNIT=85,STATUS='NEW',DEFAULTFILE='.OUT',FILE=FILEX)
C      IKEY=LIB$SPAWN('$ SET TERM/WIDTH=132')
C      WRITE(85,5001)
C      write(6,50041) filex
C      WRITE(85,50041)filex
C50041 format(1x,'THE FOLLOWING FILE IS BEING USED 'a)
C      WRITE(6,54099)
      WRITE(85,54099)
54099 FORMAT(' ****** NOTE FROM EXPERIENCE, UNDER THE FIGURE OF',
     11X,'MERIT, LOOK FOR HIGH VALUES--THESE COULD BE CORRECT'//
     110X,'VAX VERSION SEPT 24, 1987--LOUER"S VERSION SEPT 30, 1982'//
     110x,'A FOM > 25 SHOULD PROBABLY BE GIVEN FURTHER THOUGHT.'/
     210X,'WITH THESE SOLUTIONS, IT IS SUGGESTED THAT THE RESULTS'
     3,1X,'BE USED WITH D-GEN FOR FURTHER INVESTIGATION.'/
     110X,'I HAVE LIMITED THE NUMBER OF ANSWERS TO 10'/
     110X,'YOU MAY WANT TO GO THROUGH "TRAZER" TO SEE EQUIVALENT'
     1,1X,'CELLS (CELL REDUCTION).'//
     210X,'NOTE THAT WITH THIS CODE, ALL LINES ARE INDEXED IN'
     1  ,1X,'A PRINTED SOLUTION.'//
     41X,'****** NOTE A SUMMARY OF THE RESULTS ARE IN XXXXXXX.OUT'
     5,1X,'******'////) 
C      IKEY=LIB$WAIT(5.0) 
CD      IKEY=LIB$SPAWN('$ WAIT 00:00:05')
C     L E C T U R E   D E S   D O N N E E S                                     
      READ(5,525,ERR=9977,END=9977)(IDENT(K),K=1,16)                     
      READ(5,600,ERR=9977,END=9977)N,ITYPE,JC,JT,JH,JO,JM                    
      READ(5,601,ERR=9977,END=9977)AMAX,BMAX,CMAX,VOLMIN,VOLMAX,BEMIN,
     1     BEMAX                       
      READ(5,602,ERR=9977,END=9977)WAVE,POIMOL,DENS,DELDEN                   
      READ(5,603,ERR=9977,END=9977)(D(I),EPSIL(I),I=1,N)                  
C                                                                               
C     P R E P A R A T I O N   E T   E C R I T U R E                             
C                 D E S   D O N N E E S                                         
C                                                                               
      WAVE2=WAVE/2.                                                             
      PAS=0.5                                                                   
      COEF4=1.E-04                                                              
      COEFF=1.E+08                                                              
      PI=4.*ATAN(1.)                                                            
      PIRAD=PI/180.                                                             
      PIDEG=180./PI                                                             
      WRITE(6,515)(IDENT(K),K=1,16)                                             
      WRITE(85,515) (IDENT(K),K=1,16)
      IF(ITYPE.EQ.1) THEN
                     WRITE(6,11515)
                     WRITE(85,11515)
11515 FORMAT('                      DATA IS IN THETA')
      ENDIF
      IF(ITYPE.EQ.2) THEN
                     WRITE(6,11516)
                     WRITE(85,11516)
11516 FORMAT('                      DATA IS IN 2 THETA')
      ENDIF
      IF(ITYPE.EQ.3) THEN
                     WRITE(6,11517)
                     WRITE(85,11517)
11517 FORMAT('                      DATA IS IN D')    
      ENDIF
      IF(ITYPE.EQ.4) THEN
                     WRITE(6,11518)
                     WRITE(85,11518)
11518 FORMAT('                      DATA IS IN Q')
      ENDIF
      WRITE(6,506)                                                              
      DTH=0.                                                                    
      GOTO(110,112,114,116),ITYPE                                               
  116 WRITE(6,501)                                                              
      WRITE(6,511)(D(I),EPSIL(I),I=1,N)                                         
      WRITE(85,511)(D(I),EPSIL(I),I=1,N)
      DO 118 I=1,N                                                              
      Q(I)=D(I)*COEF4                                                           
      EPSQ(I)=EPSIL(I)*COEF4                                                    
      D(I)=1./SQRT(Q(I))                                                        
      QTERME=EPSQ(I)/(2.*Q(I))                                                  
      EPSIL(I)=D(I)*QTERME                                                      
      ANGULO=ASIN(WAVE2/D(I))                                                   
      TH(I)=2.*PIDEG*ANGULO                                                     
      DTH=DTH+(QTERME*PIDEG*TAN(ANGULO))                                        
  118 CONTINUE                                                                  
      GOTO 122                                                                  
  114 WRITE(6,502)                                                              
      WRITE(6,510)(D(I),EPSIL(I),I=1,N)                                         
      WRITE(85,510)(D(I),EPSIL(I),I=1,N)                                      
      DO 124 I=1,N                                                              
      ANGULO=ASIN(WAVE2/D(I))                                                   
      TH(I)=2.*PIDEG*ANGULO                                                     
      DTH=DTH+(2.*TAN(ANGULO)*EPSIL(I)*PIDEG/D(I))                              
  124 CONTINUE                                                                  
      GOTO 126                                                                  
  112 WRITE(6,503)                                                              
      WRITE(6,511)(D(I),EPSIL(I),I=1,N)                                         
      WRITE(85,511)(D(I),EPSIL(I),I=1,N)
      DOS=1.                                                                    
      DOSINV=.5                                                                 
      GOTO 128                                                                  
  110 WRITE(6,504)                                                              
      WRITE(6,511)(D(I),EPSIL(I),I=1,N)                                         
      WRITE(85,511)(D(I),EPSIL(I),I=1,N)
      DOS=2.                                                                    
      DOSINV=1.                                                                 
  128 DO 132 I=1,N                                                              
      TH(I)=DOS*D(I)                                                            
      ANGULO=PIRAD*DOSINV*D(I)                                                  
      D(I)=WAVE2/SIN(ANGULO)                                                    
      DTH=DTH+DOS*EPSIL(I)                                                      
      EPSIL(I)=D(I)*EPSIL(I)*PIRAD*DOSINV/TAN(ANGULO)                           
  132 CONTINUE                                                                  
  126 DO 134 I=1,N                                                              
      PR=D(I)*D(I)                                                              
      Q(I)=1./PR                                                                
      EPSQ(I)=2.*EPSIL(I)/(PR*D(I))                                             
  134 CONTINUE                                                                  
  122 DTH=DTH/N                                                                 
      DO 101 I=1,N                                                              
      KQ(I)=COEFF*Q(I)                                                          
      KEPSQ(I)=COEFF*EPSQ(I)                                                    
  101 CONTINUE                                                                  
      IF (JM.EQ.0) GOTO 201                                                     
      IF (BEMIN.EQ.0.0) BEMIN=90.                                               
      IF (BEMAX.EQ.0.0) BEMAX=130.                                              
      WRITE(6,524)AMAX,BEMIN,VOLMIN,BMAX,BEMAX,VOLMAX,CMAX                      
      WRITE(85,524)AMAX,BEMIN,VOLMIN,BMAX,BEMAX,VOLMAX,CMAX  
      GOTO 202                                                                  
  201 WRITE(6,526)AMAX,VOLMIN,BMAX,VOLMAX,CMAX                                  
      WRITE(85,526)AMAX,VOLMIN,BMAX,VOLMAX,CMAX                  
  202 WRITE(6,527)WAVE                                                          
C                                                                               
C     A V E C   O U   S A N S   D E N S I T E ?                                 
C                                                                               
      KDENS=1                                                                   
      X=POIMOL*DENS*DELDEN                                                      
      IF(X.LE.1.0E-04)GO TO 70                                                  
      KDENS=2                                                                   
      WRITE(6,528)POIMOL,DENS,DELDEN                                            
   70 WRITE(6,22)                                                               
C                                                                               
C     F O R M A T S                                                             
C                                                                               
    4 FORMAT(/1H0,10X,18HPOIDS MOLECULAIRE=,F9.4,10X,8HDENSITE=,F7.4,10X        
     1,10HPRECISION=,F7.4)                                                      
  515 FORMAT(1H1/1H0,26X,68(1H-)/1H ,26X,2H! ,16A4,2H !/1H ,26X,68(1H-))        
  506 FORMAT(/45X,30HDONNEES FOURNIES A L EXECUTION/45X,30(1H*))                
    7 FORMAT(//1H0,55X,10(1H*)////1H0,10X,11HA MAXIMUM =,F6.2/1H0,10X,11        
     1HB MAXIMUM =,F6.2/1H0,10X,11HC MAXIMUM =,F6.2/1H0,10X,16HVOLUME MI        
     2NIMUM =F9.2/1H0,10X,16HVOLUME MAXIMUM =,F9.2)                             
    8 FORMAT(5X,'SYSTEME CUBIQUE'/5X,15('-'))                                   
    9 FORMAT(52X,'PAS  D INDEXATION')                                           
   10 FORMAT(5X,'SYSTEME TETRAGONAL'/5X,18('-'))
   11 FORMAT(5X,'SYSTEME HEXAGONAL'/5X,17('-'))                              
   12 FORMAT(5X,'SYSTEME ORTHORHOMBIQUE'/5X,22('-'))                            
   13 FORMAT(1H0,5X,83HFIN DE LA RECHERCHE D UNE SOLUTION CUBIQUE,TETRAG        
     1ONALE,HEXAGONALE OU ORTHORHOMBIQUE/1H ,5X,83(1H*))                        
  510 FORMAT(1H ,33X,F10.4,36X,F10.4)                                           
   16 FORMAT(5X,'SYSTEME MONOCLINIQUE'/5X,20('-'))                              
   17 FORMAT(//10X,'DOMAINE VOLUME EXAMINE (TOLERANCE SUPERIEURE)'/1H           
     1,10X,45('*')//1H0,20X,'VOLUME MINIMUM=',F9.2,'A**3',10X,                  
     2'VOLUME MAXIMUN=',F9.2,'A**3'//)                                          
   22 FORMAT(///30X,60('*')/30X,'**',56X,'**'/30X,'**  ATTENTION  : VOS         
     SDONNEES SONT-ELLES IRREPROCHABLES ?  **'/30X,'**',56X,'**'/30X,           
     S60('*'))                                                                  
   23 FORMAT(1H0,5X,47HFIN DE LA RECHERCHE D UNE SOLUTION MONOCLINIQUE/1        
     1H ,5X,47(1H*))                                                            
   24 FORMAT(/1H0,10X,14HBETA MINIMUM =,F9.2/1H0,10X,14HBETA MAXIMUM =,F        
     19.2)                                                                      
  525 FORMAT(16A4)                                                              
  600 FORMAT(7I2)                                                               
  601 FORMAT(7F10.2)                                                            
  602 FORMAT(4F10.5)                                                            
  603 FORMAT(2F10.4)                                                            
  501 FORMAT(35X,'Q-UNITES',41X,'ERREUR'/31X,'EXPERIMENTALES',                  
     S       34X,'EXPERIMENTALE'/)                                              
  502 FORMAT(28X,'DISTANCE INTERRETICULAIRE',31X,'ERREUR'/                      
     S       4X,2(31X,'EXPERIMENTALE')/)                                        
  503 FORMAT(38X,'2-THETA',39X,'ERREUR'/35X,'EXPERIMENTAL',                     
     S       33X,'EXPERIMENTALE'/)                                              
  504 FORMAT(39X,'THETA',40X,'ERREUR'/35X,'EXPERIMENTAL',                       
     S       33X,'EXPERIMENTALE'/)                                              
  511 FORMAT(33X,F10.3,36X,F10.3)                                               
  524 FORMAT(//55X,10('*')//15X,9('-'),' LIMITES SUR LES PARAMETRES ',          
     S       10('-'),5x,5('-'),' LIMITES SUR LE VOLUME ',6('-')/                
     S       15X,'!',45X,'!',5X,'!',32X,'!'/15X,'!',4X,'A MAXIMUM =',           
     S       F6.2,24X,'!',5X,'!',32X,'!'/15X,'!',19X,'BETA MINIMUM =',          
     S       F8.2,4X,'!',5X,'!',4X,'VOLUME MINIMUM =',F8.2,4X,'!'/              
     S       15X,'!',4X,'B MAXIMUM =',F6.2,24X,'!',5X,'!',32X,'!'/              
     S       15X,'!',19X,'BETA MAXIMUM =',F8.2,4X,'!',5X,'!',4X,                
     S       'VOLUME MAXIMUM =',F8.2,4X,'!'/15X,'!',4X,'C MAXIMUM =',           
     S       F6.2,24X,'!',5X,'!',32X,'!'/15X,'!',45X,'!',5X,'!',32X,            
     S       '!'/15X,47('-'),5X,34('-')//)                                      
  526 FORMAT(//55X,10('*')//24X,'-- LIMITES SUR LES PARAMETRES --',             
     S       8X,6('-'),' LIMITES SUR LE VOLUME ',5('-')/                        
     S       24X,'!',30X,'!',8X,'!',32X,'!'/24X,'!',7X,'A MAXIMUM =',           
     S       F6.2,6X,'!',8X,'!',32X,'!'/24X,'!',30X,'!',8X,'!',4X,              
     S       'VOLUME MINIMUM =',F8.2,4X,'!'/24X,'!',7X,'B MAXIMUM =',           
     S       F6.2,6X,'!',8X,'!',32X,'!'/24X,'!',30X,'!',8X,'!',4X,              
     S       'VOLUME MAXIMUM =',F8.2,4X,'!'/24X,'!',7X,'C MAXIMUM =',           
     S       F6.2,6X,'!',8X,'!',32X,'!'/24X,'!',30X,'!',8X,'!',32X,             
     S       '!'/24X,32('-'),8X,34('-')//)                                      
  527 FORMAT(50X,'LONGUEUR  D ONDE =',F9.6)                                     
  528 FORMAT(//5X,'POIDS MOLECULAIRE =',F9.4,5X,                                
     S       'DENSITE OBSERVEE =',F8.4,' G.CM(-3)',5X,                          
     S       'ERREUR S/DENSITE =',F8.4,' G.CM(-3)'/)                            
C                                                                               
C     P A R T I E   C O M M U N E                                               
C                                                                               
      V=VOLMAX                                                                  
      VV=V                                                                      
      PAS2=PAS/2.                                                               
      PAS4=PAS/4.                                                               
      PAS8=PAS/8.                                                               
      PAS16=PAS/16.                                                             
      PAS32=PAS/32.                                                             
      PAS64=PAS/64.                                                             
      VSUPM=VOLMIN                                                              
      AMIN=D(1)-EPSIL(1)                                                        
      ZB=2.*EPSIL(2)+EPSIL(1)                                                   
      ZZ=ABS(2.*D(2)-D(1))                                                      
      IF(ZZ-ZB)19,19,18                                                         
   19 BMIN=D(3)-EPSIL(3)                                                        
      GO TO 27                                                                  
   18 BMIN=D(2)-EPSIL(2)                                                        
   27 CMIN=2.5                                                                  
      NA=(AMAX-AMIN)/PAS+1.                                                     
      NB=(BMAX-BMIN)/PAS+1.                                                     
      NC=(CMAX-CMIN)/PAS+1.                                                     
      IF(KDENS.EQ.2)GO TO 3001                                                  
C                                                                               
C     T R A I T E M E N T   S A N S   D E N S I T E                             
C                                                                               
      VSUP=VOLMIN                                                               
      PASVOL=400.                                                               
      NVOL=(VOLMAX-VOLMIN)/PASVOL                                               
      NVOL=NVOL+1                                                               
      IF(JC.EQ.0.AND.JT.EQ.0.AND.JH.EQ.0.AND.JO.EQ.0)GO TO 2100                 
      DO 1000 KVOL=1,NVOL                                                       
      V1=VV                                                                     
      VINF=VSUP                                                                 
      VSUP=VINF+PASVOL                                                          
      IF(VSUP.GT.VOLMAX)VSUP=VOLMAX                                             
      IF(VINF.GE.VOLMAX)GO TO 1100                                              
      WRITE(6,17)VINF,VSUP
      WRITE(85,17)VINF,VSUP
      IF(JC.EQ.0)GO TO 21                                                       
      WRITE(6,8)                                                                
      CALL CUBIQU(AMIN,AMAX,NA)                                                 
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)20,20,21                                                   
   20 WRITE(6,9)                                                                
   21 IF(JT.EQ.0)GO TO 31                                                       
      VV=V                                                                      
      WRITE(6,10)                                                               
      ICHOIX=-1                                                                 
      CALL TETHEX(CMIN,BMIN,NC,ICHOIX)                                          
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)30,30,31                                                   
   30 WRITE(6,9)                                                                
   31 IF(JH.EQ.0)GO TO 41                                                       
      VV=V                                                                      
      WRITE(6,11)                                                               
      ICHOIX=1                                                                  
      CALL TETHEX(CMIN,BMIN,NC,ICHOIX)                                          
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)40,40,41                                                   
   40 WRITE(6,9)                                                                
   41 IF(JO.EQ.0)GO TO 51                                                       
      VV=V                                                                      
      WRITE(6,12)                                                               
      CALL ORTHOR(AMIN,BMIN,CMIN,NA,NB,NC)                                      
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)50,50,51                                                   
   50 WRITE(6,9)                                                                
   51 DIFVOL=ABS(V-V1)                                                          
      IF(DIFVOL-0.02)1000,1000,1100                                             
 1000 CONTINUE                                                                  
 1100 WRITE(6,13)                                                               
      IF(JM.EQ.0)GO TO 2200                                                     
      VV=V                                                                      
 2100 NA=(AMAX-BMIN)/PAS+1.                                                     
      NB=(BMAX-CMIN)/PAS+1.                                                     
      BEMIN=BEMIN*PIRAD                                                         
      BEMAX=BEMAX*PIRAD                                                         
      VSUP=VSUPM                                                                
      DO 2000 KVOL=1,NVOL                                                       
      V1=VV                                                                     
      VINF=VSUP                                                                 
      VSUP=VINF+PASVOL                                                          
      IF(VSUP.GT.VOLMAX)VSUP=VOLMAX                                             
      IF(VINF.GT.VOLMAX)GO TO 2300                                              
      WRITE(6,17)VINF,VSUP                                                      
      WRITE(85,17)VINF,VSUP
      WRITE(6,16)                                                               
      CALL MONOC(BMIN,NA,NB,CMIN,NC,BEMAX,BEMIN)                                
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)60,60,61                                                   
   60 WRITE(6,9)                                                                
   61 DIFVOL=ABS(V-V1)                                                          
      IF(DIFVOL-0.02)2000,2000,2300                                             
 2000 CONTINUE                                                                  
 2300 WRITE(6,23)                                                               
      GO TO 2200                                                                
C                                                                               
C     T R A I T E M E N T   A V E C   D E N S I T E                             
C                                                                               
 3001 AVOG=0.602253                                                             
      RAP=POIMOL/AVOG                                                           
      VUNITP=RAP/(DENS-DELDEN)                                                  
      VUNITM=RAP/(DENS+DELDEN)                                                  
      IF(JC.EQ.0.AND.JT.EQ.0.AND.JH.EQ.0.AND.JO.EQ.0)GO TO 3200                 
      DO 3000 KZ=1,50                                                           
      V1=VV                                                                     
      VINF=KZ*VUNITM                                                            
      VSUP=KZ*VUNITP                                                            
      IF (VSUP.LT.VOLMIN)GO TO 3000                                             
      IF(VINF.GE.VOLMAX)GO TO 3100                                              
      WRITE(6,17)VINF,VSUP                                                      
      WRITE(85,17)VINF,VSUP
      IF(JC.EQ.0)GO TO 121                                                      
      WRITE(6,8)                                                                
      CALL CUBIQU(AMIN,AMAX,NA)                                                 
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)120,120,121                                                
  120 WRITE(6,9)                                                                
  121 IF(JT.EQ.0)GO TO 131                                                      
      VV=V                                                                      
      WRITE(6,10)                                                               
      ICHOIX=-1                                                                 
      CALL TETHEX(CMIN,BMIN,NC,ICHOIX)                                          
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)130,130,131                                                
  130 WRITE(6,9)                                                                
  131 IF(JH.EQ.0)GO TO 141                                                      
      VV=V                                                                      
      WRITE(6,11)                                                               
      ICHOIX=1                                                                  
      CALL TETHEX(CMIN,BMIN,NC,ICHOIX)                                          
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)140,140,141                                                
  140 WRITE(6,9)                                                                
  141 IF(JO.EQ.0)GO TO 3000                                                     
      VV=V                                                                      
      WRITE(6,12)                                                               
      CALL ORTHOR(AMIN,BMIN,CMIN,NA,NB,NC)                                      
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)150,150,151                                                
  150 WRITE(6,9)                                                                
  151 DIFVOL=ABS(V-V1)                                                          
      IF(DIFVOL-0.02)3000,3000,3100                                             
 3000 CONTINUE                                                                  
 3100 WRITE(6,13)                                                               
      IF(JM.EQ.0)GO TO 2200                                                     
      VV=V                                                                      
 3200 NA=(AMAX-BMIN)/PAS+1.                                                     
      NB=(BMAX-CMIN)/PAS+1.                                                     
      BEMIN=BEMIN*PIRAD                                                         
      BEMAX=BEMAX*PIRAD                                                         
      VSUP=VSUPM                                                                
      DO 4000 KZ=1,50                                                           
      V1=VV                                                                     
      VINF=KZ*VUNITM                                                            
      VSUP=KZ*VUNITP                                                            
      IF(VSUP.LT.VOLMIN)GO TO 4000                                              
      IF(VINF.GE.VOLMAX)GO TO 4300                                              
      WRITE(6,17)VINF,VSUP                                                      
      WRITE(85,17)VINF,VSUP
      WRITE(6,16)                                                               
      CALL MONOC(BMIN,NA,NB,CMIN,NC,BEMAX,BEMIN)                                
      DIFVOL=ABS(VV-V)                                                          
      IF(DIFVOL-0.02)160,160,161                                                
  160 WRITE(6,9)                                                                
  161 DIFVOL=ABS(V-V1)                                                          
      IF(DIFVOL-0.02)4000,4000,4300                                             
 4000 CONTINUE                                                                  
 4300 WRITE(6,23)                                                               
      CLOSE(UNIT=5)
      CLOSE(UNIT=85)
 2200 STOP 'END OF LOUER"S DICVOL PROGRAM'
 9977 STOP 'END OF FILE ON READ STATEMENTS'                                 
      END                                                                       
C     ------------------------------------------------------------------        
C     !                         C U B I Q U E                          !        
C     ------------------------------------------------------------------        
      SUBROUTINE CUBIQU(AMIN,AMAX,NA)                                           
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      COMMON QCAL(50,11),XH(50,11),XK(50,11),XL(50,11),AA,BB,CC,BETA            
      DIMENSION NDICH(7)                                                        
    1 FORMAT(1H0,10X,2HA=,F8.5,10X,33HVOLUME DE LA MAILLE ELEMENTAIRE =,        
     1F10.3)                                                                    
    2 FORMAT(1H ,10X,F8.5,4X,3(10X,F3.0),9X,F8.5,11X,F8.5,9X,F8.5)              
    3 FORMAT(1H0,50X,9(1H-)/)                                                   
    4 FORMAT(1H0,10X,2HA=,F8.5,10X,33HVOLUME DE LA MAILLE ELEMENTAIRE =,        
     1F10.3,5X,2HZ=,I2)                                                         
    5 FORMAT(1H0,7X,14HD EXPERIMENTAL,12X,1HH,12X,1HK,12X,1HL,10X,9HD CA        
     1LCULE,12X,5HDE-DC,12X,7HEPSILON/1H0)                                      
      DO 21 I=1,7                                                               
   21 NDICH(I)=0                                                                
      AA=1000.                                                                  
      DO 1000 IA=1,NA                                                           
      AW=IA-1                                                                   
      AMOI1=AMIN+(AW*PAS)                                                       
      IF(AMOI1.GT.AMAX)GO TO 1000                                               
      APLU1=AMOI1+PAS                                                           
      VTESTP=APLU1**3                                                           
      VTESTM=AMOI1**3                                                           
      IF(VTESTP.LT.VINF.OR.VTESTM.GT.VSUP)GO TO 1000                            
      IF(V-VTESTM)1000,20,20                                                    
   20 MH=APLU1/D(5)                                                             
      MH2=APLU1/D(N)                                                            
      KRAM1=COEFF/(AMOI1*AMOI1)                                                 
      KRAP1=COEFF/(APLU1*APLU1)                                                 
      CALL CUDHKL(KRAM1,KRAP1)                                                  
      DO 120 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 1000                                                 
  120 CONTINUE                                                                  
      NDICH(1)=NDICH(1)+1                                                       
      DO 2000 IA2=1,2                                                           
      AW=IA2-1                                                                  
      AMOI2=AMOI1+(AW*PAS2)                                                     
      IF(ABS(AMOI2-AA).LT.0.25)GO TO 2000                                       
      APLU2=AMOI2+PAS2                                                          
      VTESTM=AMOI2*AMOI2*AMOI2                                                  
      IF(VSUP.LT.VTESTM)GO TO 1000                                              
      VTESTP=APLU2*APLU2*APLU2                                                  
      IF(VINF.GT.VTESTP)GO TO 2000                                              
      KRAM2=COEFF/(AMOI2*AMOI2)                                                 
      KRAP2=COEFF/(APLU2*APLU2)                                                 
      CALL CUDHKL(KRAM2,KRAP2)                                                  
      DO 121 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 2000                                                 
  121 CONTINUE                                                                  
      NDICH(2)=NDICH(2)+1                                                       
      DO 3000 IA3=1,2                                                           
      AW=IA3-1                                                                  
      AMOI3=AMOI2+(AW*PAS4)                                                     
      IF(ABS(AMOI3-AA).LT.0.25)GO TO 3000                                       
      APLU3=AMOI3+PAS4                                                          
      VTESTM=AMOI3*AMOI3*AMOI3                                                  
      IF(VSUP.LT.VTESTM)GO TO 2000                                              
      VTESTP=APLU3*APLU3*APLU3                                                  
      IF(VINF.GT.VTESTP)GO TO 3000                                              
      KRAM3=COEFF/(AMOI3*AMOI3)                                                 
      KRAP3=COEFF/(APLU3*APLU3)                                                 
      CALL CUDHKL(KRAM3,KRAP3)                                                  
      DO 122 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 3000                                                 
  122 CONTINUE                                                                  
      NDICH(3)=NDICH(3)+1                                                       
      DO 4000 IA4=1,2                                                           
      AW=IA4-1                                                                  
      AMOI4=AMOI3+(AW*PAS8)                                                     
      IF(ABS(AMOI4-AA).LT.0.25)GO TO 4000                                       
      APLU4=AMOI4+PAS8                                                          
      VTESTM=AMOI4*AMOI4*AMOI4                                                  
      IF(VSUP.LT.VTESTM)GO TO 3000                                              
      VTESTP=APLU4*APLU4*APLU4                                                  
      IF(VINF.GT.VTESTP)GO TO 4000                                              
      KRAM4=COEFF/(AMOI4*AMOI4)                                                 
      KRAP4=COEFF/(APLU4*APLU4)                                                 
      CALL CUDHKL(KRAM4,KRAP4)                                                  
      DO 123 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 4000                                                 
  123 CONTINUE                                                                  
      NDICH(4)=NDICH(4)+1                                                       
      DO 5000 IA5=1,2                                                           
      AW=IA5-1                                                                  
      AMOI5=AMOI4+(AW*PAS16)                                                    
      IF(ABS(AMOI5-AA).LT.0.25)GO TO 5000                                       
      APLU5=AMOI5+PAS16                                                         
      VTESTM=AMOI5*AMOI5*AMOI5                                                  
      IF(VSUP.LT.VTESTM)GO TO 4000                                              
      VTESTP=APLU5*APLU5*APLU5                                                  
      IF(VINF.GT.VTESTP)GO TO 5000                                              
      KRAM5=COEFF/(AMOI5*AMOI5)                                                 
      KRAP5=COEFF/(APLU5*APLU5)                                                 
      CALL CUDHKL(KRAM5,KRAP5)                                                  
      DO 124 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 5000                                                 
  124 CONTINUE                                                                  
      NDICH(5)=NDICH(5)+1                                                       
      DO 6000 IA6=1,2                                                           
      AW=IA6-1                                                                  
      AMOI6=AMOI5+(AW*PAS32)                                                    
      IF(ABS(AMOI6-AA).LT.0.25)GO TO 6000                                       
      APLU6=AMOI6+PAS32                                                         
      VTESTM=AMOI6*AMOI6*AMOI6                                                  
      IF(VSUP.LT.VTESTM)GO TO 5000                                              
      VTESTP=APLU6*APLU6*APLU6                                                  
      IF(VINF.GT.VTESTP)GO TO 6000                                              
      KRAM6=COEFF/(AMOI6*AMOI6)                                                 
      KRAP6=COEFF/(APLU6*APLU6)                                                 
      CALL CUDHKL(KRAM6,KRAP6)                                                  
      DO 125 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 6000                                                 
  125 CONTINUE                                                                  
      NDICH(6)=NDICH(6)+1                                                       
      DO 7000 IA7=1,2                                                           
      AW=IA7-1                                                                  
      AMOI7=AMOI6+(AW*PAS64)                                                    
      IF(ABS(AMOI7-AA).LT.0.25)GO TO 7000                                       
      APLU7=AMOI7+PAS64                                                         
      VTESTM=AMOI7*AMOI7*AMOI7                                                  
      IF(VSUP.LT.VTESTM)GO TO 6000                                              
      VTESTP=APLU7*APLU7*APLU7                                                  
      IF(VINF.GT.VTESTP)GO TO 7000                                              
      RAM7=1./(AMOI7*AMOI7)                                                     
      RAP7=1./(APLU7*APLU7)                                                     
      KRAM7=COEFF*RAM7                                                          
      KRAP7=COEFF*RAP7                                                          
      CALL CUDHKL(KRAM7,KRAP7)                                                  
      DO 126 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 7000                                                 
  126 CONTINUE                                                                  
      NDICH(7)=NDICH(7)+1                                                       
      AA=AMOI7+(PAS64/2.)                                                       
      BB=AA                                                                     
      CC=AA                                                                     
      BETA=90.                                                                  
      NPOSS=0                                                                   
      DO 15 I=1,N                                                               
      DO 16 J=1,6                                                               
   16 QCAL(I,J)=-1.                                                             
   15 CONTINUE                                                                  
      DO 80 M=1,MH2                                                             
      MCARH=M*M                                                                 
      MM=M+1                                                                    
      DO 90 K=1,MM                                                              
      MXK=K-1                                                                   
      MCARK=MCARH+MXK*MXK                                                       
      DO 100 L=1,K                                                              
      MXL=L-1                                                                   
      CARL=MCARK+MXL*MXL                                                        
      QQP=CARL*RAM7                                                             
      QQM=CARL*RAP7                                                             
      IF(QQM.GT.(Q(N)+EPSQ(N)))GO TO 90                                         
      NPOSS=NPOSS+1                                                             
      DO 110 I=1,N                                                              
      IF(QQP.LT.(Q(I)-EPSQ(I)))GO TO 110                                        
      IF(QQM.GT.(Q(I)+EPSQ(I)))GO TO 110                                        
      J=1                                                                       
  111 IF(QCAL(I,J))112,113,113                                                  
  113 J=J+1                                                                     
      IF(J-5)111,111,110                                                        
  112 QCAL(I,J)=CARL/AA**2                                                      
      XH(I,J)=M                                                                 
      XK(I,J)=MXK                                                               
      XL(I,J)=MXL                                                               
  110 CONTINUE                                                                  
  100 CONTINUE                                                                  
   90 CONTINUE                                                                  
   80 CONTINUE                                                                  
      DO 131 I=1,N                                                              
      IF(QCAL(I,1))7000,131,131                                                 
  131 CONTINUE                                                                  
      IND=1                                                                     
      NRIND=1                                                                   
      CALL AFFPAR(IND,NRIND,NPOSS)                                              
      WRITE(6,3)                                                                
 7000 CONTINUE                                                                  
 6000 CONTINUE                                                                  
 5000 CONTINUE                                                                  
 4000 CONTINUE                                                                  
 3000 CONTINUE                                                                  
 2000 CONTINUE                                                                  
 1000 CONTINUE                                                                  
      WRITE(6,10)(NDICH(I),I=1,7)                                               
   10 FORMAT(1H0,10X,42HNOMBRE D ITERATIONS DANS LES 7 DICHOTOMIES/1H ,1        
     10X,7(3X,I10))                                                             
      RETURN                                                                    
      END                                                                       
      SUBROUTINE CUDHKL(KRAM,KRAP)                                              
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      DO 15 I=1,N                                                               
      KQT(I)=-1                                                                 
   15 CONTINUE                                                                  
      DO 10 M=1,MH                                                              
      MCARH=M*M                                                                 
      MM=M+1                                                                    
      DO 20 K=1,MM                                                              
      MCARK=(K-1)*(K-1)+MCARH                                                   
      DO 30 L=1,K                                                               
      MCARL=(L-1)*(L-1)+MCARK                                                   
      KQQM=MCARL*KRAP                                                           
      IF((KQQM-KEPSQ(5)).GT.KQ(5))GO TO 20                                      
      KQQP=MCARL*KRAM                                                           
      IF((KQQP+KEPSQ(1)).LT.KQ(1))GO TO 30                                      
      DO 40 I=1,5                                                               
      IF((KQQM-KEPSQ(I)).GT.KQ(I))GO TO 40                                      
      IF((KQQP+KEPSQ(I)).GE.KQ(I))KQT(I)=1                                      
   40 CONTINUE                                                                  
   30 CONTINUE                                                                  
   20 CONTINUE                                                                  
   10 CONTINUE                                                                  
      DO 70 I=1,5                                                               
      IF(KQT(I).LT.0)GO TO 3000                                                 
   70 CONTINUE                                                                  
      DO 80 M=1,MH2                                                             
      MCARH=M*M                                                                 
      MM=M+1                                                                    
      DO 90 K=1,MM                                                              
      MCARK=(K-1)*(K-1)+MCARH                                                   
      DO 100 L=1,K                                                              
      MCARL=(L-1)*(L-1)+MCARK                                                   
      KQQP=MCARL*KRAM                                                           
      IF((KQQP+KEPSQ(6)).LT.KQ(6))GO TO 100                                     
      KQQM=MCARL*KRAP                                                           
      IF((KQQM-KEPSQ(N)).GT.KQ(N))GO TO 90                                      
      DO 110 I=6,N                                                              
      IF((KQQM-KEPSQ(I)).GT.KQ(I))GO TO 110                                     
      IF((KQQP+KEPSQ(I)).GE.KQ(I))KQT(I)=1                                      
  110 CONTINUE                                                                  
  100 CONTINUE                                                                  
   90 CONTINUE                                                                  
   80 CONTINUE                                                                  
 3000 RETURN                                                                    
      END                                                                       
C     ------------------------------------------------------------------        
C     !          T E T R A G O N A L   E T   H E X A G O N A L         !        
C     ------------------------------------------------------------------        
      SUBROUTINE TETHEX(CMIN,BMIN,NC,ICHOIX)                                    
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      COMMON QCAL(50,11),XH(50,11),XK(50,11),XL(50,11),AA,BB,CC,BETA            
      DIMENSION NDICH(7)                                                        
    1 FORMAT(1H0,10X,2HA=,F8.5,10X,2HC=,F8.5,10X,33HVOLUME DE LA MAILLE         
     1ELEMENTAIRE =,F10.3)                                                      
    2 FORMAT(1H ,10X,F8.5,4X,3(10X,F3.0),9X,F8.5,11X,F8.5,9X,F8.5)              
    3 FORMAT(1H0,50X,9(1H-)/)                                                   
    4 FORMAT(1H0,10X,2HA=,F8.5,10X,2HC=,F8.5,10X,33HVOLUME DE LA MAILLE         
     1ELEMENTAIRE =,F10.3,5X,2HZ=,I2)                                           
    5 FORMAT(1H0,7X,14HD EXPERIMENTAL,12X,1HH,12X,1HK,12X,1HL,10X,9HD CA        
     1LCULE,12X,5HDE-DC,12X,7HEPSILON/1H0)                                      
      DO 20 I=1,7                                                               
   20 NDICH(I)=0                                                                
      AA=1000.                                                                  
      CC=1000.                                                                  
      RA=4./3.                                                                  
      CONST=COEFF                                                               
      CONTR=D(1)-EPSIL(1)                                                       
       CONTR1=CONTR                                                             
      COVOL=1.                                                                  
      IF(ICHOIX)79,79,78                                                        
   78 CONST=RA*COEFF                                                            
      COVOL=0.8660254                                                           
      CONTR=CONTR/COVOL                                                         
   79 DO 1000 IA=1,NC                                                           
      AW=IA-1                                                                   
      AMOI1=CMIN+(AW*PAS)                                                       
      APLU1=AMOI1+PAS                                                           
      MH=(APLU1/D(5))+1.                                                        
      MH2=(APLU1/D(N))+1.                                                       
      KRAM1=CONST/(AMOI1*AMOI1)                                                 
      KRAP1=CONST/(APLU1*APLU1)                                                 
      DO 1100 IC=1,NC                                                           
      CW=IC-1                                                                   
      CMOI1=CMIN+(CW*PAS)                                                       
      CPLU1=CMOI1+PAS                                                           
      IF(APLU1.LT.CONTR.AND.CPLU1.LT.CONTR1)GO TO 1100                          
      IF(CMOI1.GT.APLU1.AND.APLU1.LT.BMIN)GO TO 1100                            
      VTESTM=AMOI1*AMOI1*CMOI1*COVOL                                            
      IF(V.LT.VTESTM)GO TO 1000                                                 
      IF(VSUP.LT.VTESTM)GO TO 1000                                              
      VTESTP=APLU1*APLU1*CPLU1*COVOL                                            
      IF(VINF.GT.VTESTP)GO TO 1100                                              
      ML=(CPLU1/D(5))+1.                                                        
      ML2=(CPLU1/D(N))+1.                                                       
      KRCM1=COEFF/(CMOI1*CMOI1)                                                 
      KRCP1=COEFF/(CPLU1*CPLU1)                                                 
      CALL TEDHKL(KRAM1,KRAP1,KRCM1,KRCP1,ICHOIX)                               
      DO 120 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 1100                                                 
  120 CONTINUE                                                                  
      NDICH(1)=NDICH(1)+1                                                       
      DO 2000 IA2=1,2                                                           
      AW=IA2-1                                                                  
      AMOI2=AMOI1+(AW*PAS2)                                                     
      APLU2=AMOI2+PAS2                                                          
      KRAM2=CONST/(AMOI2*AMOI2)                                                 
      KRAP2=CONST/(APLU2*APLU2)                                                 
      DO 2100 IC2=1,2                                                           
      CW=IC2-1                                                                  
      CMOI2=CMOI1+(CW*PAS2)                                                     
      ATEST=ABS(AMOI2-AA)                                                       
      CTEST=ABS(CMOI2-CC)                                                       
      IF(ATEST.LT.0.25.AND.CTEST.LT.0.25)GO TO 2100                             
      VTESTM=AMOI2*AMOI2*CMOI2*COVOL                                            
      IF(VSUP.LT.VTESTM)GO TO 2000                                              
      CPLU2=CMOI2+PAS2                                                          
      VTESTP=APLU2*APLU2*CPLU2*COVOL                                            
      IF(VINF.GT.VTESTP)GO TO 2100                                              
      KRCM2=COEFF/(CMOI2*CMOI2)                                                 
      KRCP2=COEFF/(CPLU2*CPLU2)                                                 
      CALL TEDHKL(KRAM2,KRAP2,KRCM2,KRCP2,ICHOIX)                               
      DO 121 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 2100                                                 
  121 CONTINUE                                                                  
      NDICH(2)=NDICH(2)+1                                                       
      DO 3000 IA3=1,2                                                           
      AW=IA3-1                                                                  
      AMOI3=AMOI2+(AW*PAS4)                                                     
      APLU3=AMOI3+PAS4                                                          
      KRAM3=CONST/(AMOI3*AMOI3)                                                 
      KRAP3=CONST/(APLU3*APLU3)                                                 
      DO 3100 IC3=1,2                                                           
      CW=IC3-1                                                                  
      CMOI3=CMOI2+(CW*PAS4)                                                     
      ATEST=ABS(AMOI3-AA)                                                       
      CTEST=ABS(CMOI3-CC)                                                       
      IF(ATEST.LT.0.25.AND.CTEST.LT.0.25)GO TO 3100                             
      VTESTM=AMOI3*AMOI3*CMOI3*COVOL                                            
      IF(VSUP.LT.VTESTM)GO TO 3000                                              
      CPLU3=CMOI3+PAS4                                                          
      VTESTP=APLU3*APLU3*CPLU3*COVOL                                            
      IF(VINF.GT.VTESTP)GO TO 3100                                              
      KRCM3=COEFF/(CMOI3*CMOI3)                                                 
      KRCP3=COEFF/(CPLU3*CPLU3)                                                 
      CALL TEDHKL(KRAM3,KRAP3,KRCM3,KRCP3,ICHOIX)                               
      DO 122 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 3100                                                 
  122 CONTINUE                                                                  
      NDICH(3)=NDICH(3)+1                                                       
      DO 4000 IA4=1,2                                                           
      AW=IA4-1                                                                  
      AMOI4=AMOI3+(AW*PAS8)                                                     
      APLU4=AMOI4+PAS8                                                          
      KRAM4=CONST/(AMOI4*AMOI4)                                                 
      KRAP4=CONST/(APLU4*APLU4)                                                 
      DO 4100 IC4=1,2                                                           
      CW=IC4-1                                                                  
      CMOI4=CMOI3+(CW*PAS8)                                                     
      ATEST=ABS(AMOI4-AA)                                                       
      CTEST=ABS(CMOI4-CC)                                                       
      IF(ATEST.LT.0.25.AND.CTEST.LT.0.25)GO TO 4100                             
      VTESTM=AMOI4*AMOI4*CMOI4*COVOL                                            
      IF(VSUP.LT.VTESTM)GO TO 4000                                              
      CPLU4=CMOI4+PAS8                                                          
      VTESTP=APLU4*APLU4*CPLU4*COVOL                                            
      IF(VINF.GT.VTESTP)GO TO 4100                                              
      KRCM4=COEFF/(CMOI4*CMOI4)                                                 
      KRCP4=COEFF/(CPLU4*CPLU4)                                                 
      CALL TEDHKL(KRAM4,KRAP4,KRCM4,KRCP4,ICHOIX)                               
      DO 123 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 4100                                                 
  123 CONTINUE                                                                  
      NDICH(4)=NDICH(4)+1                                                       
      DO 5000 IA5=1,2                                                           
      AW=IA5-1                                                                  
      AMOI5=AMOI4+(AW*PAS16)                                                    
      APLU5=AMOI5+PAS16                                                         
      KRAM5=CONST/(AMOI5*AMOI5)                                                 
      KRAP5=CONST/(APLU5*APLU5)                                                 
      DO 5100 IC5=1,2                                                           
      CW=IC5-1                                                                  
      CMOI5=CMOI4+(CW*PAS16)                                                    
      ATEST=ABS(AMOI5-AA)                                                       
      CTEST=ABS(CMOI5-CC)                                                       
      IF(ATEST.LT.0.25.AND.CTEST.LT.0.25)GO TO 5100                             
      VTESTM=AMOI5*AMOI5*CMOI5*COVOL                                            
      IF(VSUP.LT.VTESTM)GO TO 5000                                              
      CPLU5=CMOI5+PAS16                                                         
      VTESTP=APLU5*APLU5*CPLU5*COVOL                                            
      IF(VINF.GT.VTESTP)GO TO 5100                                              
      KRCM5=COEFF/(CMOI5*CMOI5)                                                 
      KRCP5=COEFF/(CPLU5*CPLU5)                                                 
      CALL TEDHKL(KRAM5,KRAP5,KRCM5,KRCP5,ICHOIX)                               
      DO 124 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 5100                                                 
  124 CONTINUE                                                                  
      NDICH(5)=NDICH(5)+1                                                       
      DO 6000 IA6=1,2                                                           
      AW=IA6-1                                                                  
      AMOI6=AMOI5+(AW*PAS32)                                                    
      APLU6=AMOI6+PAS32                                                         
      KRAM6=CONST/(AMOI6*AMOI6)                                                 
      KRAP6=CONST/(APLU6*APLU6)                                                 
      DO 6100 IC6=1,2                                                           
      CW=IC6-1                                                                  
      CMOI6=CMOI5+(CW*PAS32)                                                    
      ATEST=ABS(AMOI6-AA)                                                       
      CTEST=ABS(CMOI6-CC)                                                       
      IF(ATEST.LT.0.25.AND.CTEST.LT.0.25)GO TO 6100                             
      VTESTM=AMOI6*AMOI6*CMOI6*COVOL                                            
      IF(VSUP.LT.VTESTM)GO TO 6000                                              
      CPLU6=CMOI6+PAS32                                                         
      VTESTP=APLU6*APLU6*CPLU6*COVOL                                            
      IF(VINF.GT.VTESTP)GO TO 6100                                              
      KRCM6=COEFF/(CMOI6*CMOI6)                                                 
      KRCP6=COEFF/(CPLU6*CPLU6)                                                 
      CALL TEDHKL(KRAM6,KRAP6,KRCM6,KRCP6,ICHOIX)                               
      DO 125 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 6100                                                 
  125 CONTINUE                                                                  
      NDICH(6)=NDICH(6)+1                                                       
      DO 7000 IA7=1,2                                                           
      AW=IA7-1                                                                  
      AMOI7=AMOI6+(AW*PAS64)                                                    
      APLU7=AMOI7+PAS64                                                         
      RAM7=1./(AMOI7*AMOI7)                                                     
      RAP7=1./(APLU7*APLU7)                                                     
      KRAM7=CONST*RAM7                                                          
      KRAP7=CONST*RAP7                                                          
      DO 7100 IC7=1,2                                                           
      CW=IC7-1                                                                  
      CMOI7=CMOI6+(CW*PAS64)                                                    
      ATEST=ABS(AMOI7-AA)                                                       
      CTEST=ABS(CMOI7-CC)                                                       
      IF(ATEST.LT.0.25.AND.CTEST.LT.0.25)GO TO 7100                             
      VTESTM=AMOI7*AMOI7*CMOI7*COVOL                                            
      IF(VSUP.LT.VTESTM)GO TO 7000                                              
      CPLU7=CMOI7+PAS64                                                         
      VTESTP=APLU7*APLU7*CPLU7*COVOL                                            
      IF(VINF.GT.VTESTP)GO TO 7100                                              
      RCM7=1./(CMOI7*CMOI7)                                                     
      RCP7=1./(CPLU7*CPLU7)                                                     
      KRCM7=COEFF*RCM7                                                          
      KRCP7=COEFF*RCP7                                                          
      CALL TEDHKL(KRAM7,KRAP7,KRCM7,KRCP7,ICHOIX)                               
      DO 126 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 7100                                                 
  126 CONTINUE                                                                  
      NDICH(7)=NDICH(7)+1                                                       
      AA=AMOI7+(PAS64/2.)                                                       
      CC=CMOI7+(PAS64/2.)                                                       
      BB=AA                                                                     
      BETA=90.                                                                  
      NPOSS=0                                                                   
      DO 15 I=1,N                                                               
      DO 16 J=1,6                                                               
   16 QCAL(I,J)=-1.                                                             
   15 CONTINUE                                                                  
      DO 80 M=1,MH2                                                             
      XXH=M-1                                                                   
      CARH=XXH**2                                                               
      DO 90 K=1,M                                                               
      XXK=K-1                                                                   
      CARK=CARH+XXK**2                                                          
      DO 100 L=1,ML2                                                            
      MXL=L-1                                                                   
      CARL=MXL*MXL                                                              
      IF(ICHOIX)300,300,301                                                     
  300 QQM=(CARK*RAP7)+(CARL*RCP7)                                               
      QQP=(CARK*RAM7)+(CARL*RCM7)                                               
      IF(QQM.GT.(Q(N)+EPSQ(N)))GO TO 90                                         
      NPOSS=NPOSS+1                                                             
      DO 110 I=1,N                                                              
      IF(QQP.LT.(Q(I)-EPSQ(I)))GO TO 110                                        
      IF(QQM.GT.(Q(I)+EPSQ(I)))GO TO 110                                        
      J=1                                                                       
  111 IF(QCAL(I,J))112,113,113                                                  
  113 J=J+1                                                                     
      IF(J-5)111,111,110                                                        
  112 QCAL(I,J)=(CARK/AA**2)+(CARL/CC**2)                                       
      XH(I,J)=XXH                                                               
      XK(I,J)=XXK                                                               
      XL(I,J)=MXL                                                               
  110 CONTINUE                                                                  
      GO TO 100                                                                 
  301 SOM=RA*(CARK+XXH*XXK)                                                     
      QQM=(SOM*RAP7)+(CARL*RCP7)                                                
      QQP=(SOM*RAM7)+(CARL*RCM7)                                                
      IF(QQM.GT.(Q(N)+EPSQ(N)))GO TO 90                                         
      NPOSS=NPOSS+1                                                             
      DO 115 I=1,N                                                              
      IF(QQP.LT.(Q(I)-EPSQ(I)))GO TO 115                                        
      IF(QQM.GT.(Q(I)+EPSQ(I)))GO TO 115                                        
      J=1                                                                       
  116 IF(QCAL(I,J))117,118,118                                                  
  118 J=J+1                                                                     
      IF(J-5)116,116,115                                                        
  117 QCAL(I,J)=(SOM/(AA*AA))+(CARL/(CC*CC))                                    
      XH(I,J)=XXH                                                               
      XK(I,J)=XXK                                                               
      XL(I,J)=MXL                                                               
  115 CONTINUE                                                                  
  100 CONTINUE                                                                  
   90 CONTINUE                                                                  
   80 CONTINUE                                                                  
      DO 131 I=1,N                                                              
      IF(QCAL(I,1))7100,131,131                                                 
  131 CONTINUE                                                                  
      IND=2                                                                     
      NRIND=2                                                                   
      IF(ICHOIX.GT.0)IND=3                                                      
      CALL AFFPAR(IND,NRIND,NPOSS)                                              
      WRITE(6,3)                                                                
 7100 CONTINUE                                                                  
 7000 CONTINUE                                                                  
 6100 CONTINUE                                                                  
 6000 CONTINUE                                                                  
 5100 CONTINUE                                                                  
 5000 CONTINUE                                                                  
 4100 CONTINUE                                                                  
 4000 CONTINUE                                                                  
 3100 CONTINUE                                                                  
 3000 CONTINUE                                                                  
 2100 CONTINUE                                                                  
 2000 CONTINUE                                                                  
 1100 CONTINUE                                                                  
 1000 CONTINUE                                                                  
      WRITE(6,10)(NDICH(I),I=1,7)                                               
   10 FORMAT(1H0,10X,42HNOMBRE D ITERATIONS DANS LES 7 DICHOTOMIES/1H ,1        
     10X,7(3X,I10))                                                             
      RETURN                                                                    
      END                                                                       
      SUBROUTINE TEDHKL(KRAM,KRAP,KRCM,KRCP,ICHOIX)                             
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      DO 15 I=1,N                                                               
      KQT(I)=-1                                                                 
   15 CONTINUE                                                                  
      DO 10 M=1,MH                                                              
      MM=M-1                                                                    
      MCARH=MM*MM                                                               
      DO 20 K=1,M                                                               
      KK=K-1                                                                    
      MCARK=MCARH+KK*KK                                                         
      DO 30 L=1,ML                                                              
      MCARL=(L-1)*(L-1)                                                         
      IF(ICHOIX)200,200,201                                                     
  200 KQQM=(MCARK*KRAP)+(MCARL*KRCP)                                            
      KQQP=(MCARK*KRAM)+(MCARL*KRCM)                                            
      GO TO 202                                                                 
  201 KSOM=MCARK+(MM*KK)                                                        
      KQQM=(KSOM*KRAP)+(MCARL*KRCP)                                             
      KQQP=(KSOM*KRAM)+(MCARL*KRCM)                                             
  202 IF((KQQM-KEPSQ(5)).GT.KQ(5))GO TO 20                                      
      IF((KQQP+KEPSQ(1)).LT.KQ(1))GO TO 30                                      
      DO 40 I=1,5                                                               
      IF((KQQM-KEPSQ(I)).GT.KQ(I))GO TO 40                                      
      IF((KQQP+KEPSQ(I)).GE.KQ(I))KQT(I)=1                                      
   40 CONTINUE                                                                  
   30 CONTINUE                                                                  
   20 CONTINUE                                                                  
   10 CONTINUE                                                                  
      DO 70 I=1,5                                                               
      IF(KQT(I).LT.0)GO TO 3000                                                 
   70 CONTINUE                                                                  
      DO 80 M=1,MH2                                                             
      MM=M-1                                                                    
      MCARH=MM*MM                                                               
      DO 90 K=1,M                                                               
      KK=K-1                                                                    
      MCARK=MCARH+KK*KK                                                         
      DO 100 L=1,ML2                                                            
      MCARL=(L-1)*(L-1)                                                         
      IF(ICHOIX)300,300,301                                                     
  300 KQQM=(MCARK*KRAP)+(MCARL*KRCP)                                            
      KQQP=(MCARK*KRAM)+(MCARL*KRCM)                                            
      GO TO 302                                                                 
  301 KSOM=MCARK+(MM*KK)                                                        
      KQQM=(KSOM*KRAP)+(MCARL*KRCP)                                             
      KQQP=(KSOM*KRAM)+(MCARL*KRCM)                                             
  302 IF((KQQM-KEPSQ(N)).GT.KQ(N))GO TO 90                                      
      IF((KQQP+KEPSQ(6)).LT.KQ(6))GO TO 100                                     
      DO 110 I=6,N                                                              
      IF((KQQM-KEPSQ(I)).GT.KQ(I))GO TO 110                                     
      IF((KQQP+KEPSQ(I)).GE.KQ(I))KQT(I)=1                                      
  110 CONTINUE                                                                  
  100 CONTINUE                                                                  
   90 CONTINUE                                                                  
   80 CONTINUE                                                                  
 3000 RETURN                                                                    
      END                                                                       
C     ------------------------------------------------------------------        
C     !                   O R T H O R H O M B I Q U E                  !        
C     ------------------------------------------------------------------        
      SUBROUTINE ORTHOR(AMIN,BMIN,CMIN,NA,NB,NC)                                
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      COMMON QCAL(50,11),XH(50,11),XK(50,11),XL(50,11),AA,BB,CC,BETA            
      DIMENSION NDICH(7)                                                        
    1 FORMAT(1H0,10X,2HA=,F8.5,10X,2HB=,F8.5,10X,2HC=,F8.5,10X,33HVOLUME        
     1 DE LA MAILLE ELEMENTAIRE =,F10.3)                                        
    2 FORMAT(1H ,10X,F8.5,4X,3(10X,F3.0),9X,F8.5,11X,F8.5,9X,F8.5)              
    3 FORMAT(1H0,50X,9(1H-)/)                                                   
    4 FORMAT(1H0,10X,2HA=,F8.5,10X,2HB=,F8.5,10X,2HC=,F8.5,10X,33HVOLUME        
     1 DE LA MAILLE ELEMENTAIRE =,F10.3,5X,2HZ=,I2)                             
    5 FORMAT(1H0,7X,14HD EXPERIMENTAL,12X,1HH,12X,1HK,12X,1HL,10X,9HD CA        
     1LCULE,12X,5HDE-DC,12X,7HEPSILON/1H0)                                      
      DO 20 I=1,7                                                               
   20 NDICH(I)=0                                                                
      AA=1000.                                                                  
      BB=1000.                                                                  
      CC=1000.                                                                  
CCCC      DO 1000 IA=1,NA                                             
      IA=1
1010      AW=IA-1                                                      
      AMOI1=AMIN+(AW*PAS)                                                       
      APLU1=AMOI1+PAS                                                           
      MH=(APLU1/D(5))+1.                                                        
      MH2=(APLU1/D(N))+1.                                                       
      KRAM1=COEFF/(AMOI1*AMOI1)                                                 
      KRAP1=COEFF/(APLU1*APLU1)                                                 
CCCC      DO 2000 IB=1,NB                                                    
      IB=1
2010      BW=IB-1                                                           
      BMOI1=BMIN+(BW*PAS)                                                       
      IF(APLU1-BMOI1)1000,800,800                                               
  800 BPLU1=BMOI1+PAS                                                           
      MK=(BPLU1/D(5))+1.                                                        
      MK2=(BPLU1/D(N))+1.                                                       
      KRBM1=COEFF/(BMOI1*BMOI1)                                                 
      KRBP1=COEFF/(BPLU1*BPLU1)                                                 
CCCC      DO 2500 IC=1,NC                                                     
      IC=1
2510      CW=IC-1                                                        
      CMOI1=CMIN+(CW*PAS)                                                       
      CPLU1=CMOI1+PAS                                                           
      IF(BPLU1-CMOI1)2000,801,801                                               
  801 IF(CPLU1.LE.CMIN)GO TO 2500                                               
      VTESTM=AMOI1*BMOI1*CMOI1                                                  
      IF(V.LT.VTESTM)GO TO 2000                                                 
      IF(VSUP.LT.VTESTM)GO TO 2000                                              
      VTESTP=APLU1*BPLU1*CPLU1                                                  
      IF(VINF.GT.VTESTP)GO TO 2500                                              
      ML=(CPLU1/D(5))+1.                                                        
      ML2=(CPLU1/D(N))+1.                                                       
      KRCM1=COEFF/(CMOI1*CMOI1)                                                 
      KRCP1=COEFF/(CPLU1*CPLU1)                                                 
      CALL ORDHKL(KRAM1,KRAP1,KRBM1,KRBP1,KRCM1,KRCP1)                          
      DO 120 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 2500                                                 
  120 CONTINUE                                                                  
      NDICH(1)=NDICH(1)+1                                                       
CCCC      DO 3000 IA2=1,2                                                     
      IA2=1
3010      AW=IA2-1                                                             
      AMOI2=AMOI1+(AW*PAS2)                                                     
      APLU2=AMOI2+PAS2                                                          
      KRAM2=COEFF/(AMOI2*AMOI2)                                                 
      KRAP2=COEFF/(APLU2*APLU2)                                                 
CCCC      DO 3100 IB2=1,2                                                       
      IB2=1
3110      BW=IB2-1                                                             
      BMOI2=BMOI1+(BW*PAS2)                                                     
      BPLU2=BMOI2+PAS2                                                          
      KRBM2=COEFF/(BMOI2*BMOI2)                                                 
      KRBP2=COEFF/(BPLU2*BPLU2)                                                 
CCCC      DO 3200 IC2=1,2                                                     
      IC2=1
3210      CW=IC2-1                                                            
      CMOI2=CMOI1+(CW*PAS2)                                                     
      ATEST=ABS(AMOI2-AA)                                                       
      BTEST=ABS(BMOI2-BB)                                                       
      CTEST=ABS(CMOI2-CC)                                                       
      IF(ATEST.LT.0.25.AND.BTEST.LT.0.25.AND.CTEST.LT.0.25) GO TO 3200          
      CPLU2=CMOI2+PAS2                                                          
      VTESTM=AMOI2*BMOI2*CMOI2                                                  
      IF(VSUP.LT.VTESTM)GO TO 3100                                              
      VTESTP=APLU2*BPLU2*CPLU2                                                  
      IF(VINF.GT.VTESTP)GO TO 3200                                              
      KRCM2=COEFF/(CMOI2*CMOI2)                                                 
      KRCP2=COEFF/(CPLU2*CPLU2)                                                 
      CALL ORDHKL(KRAM2,KRAP2,KRBM2,KRBP2,KRCM2,KRCP2)                          
      DO 121 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 3200                                                 
  121 CONTINUE                                                                  
      NDICH(2)=NDICH(2)+1                                                       
      DO 4000 IA3=1,2                                                           
      AW=IA3-1                                                                  
      AMOI3=AMOI2+(AW*PAS4)                                                     
      APLU3=AMOI3+PAS4                                                          
      KRAM3=COEFF/(AMOI3*AMOI3)                                                 
      KRAP3=COEFF/(APLU3*APLU3)                                                 
      DO 4100 IB3=1,2                                                           
      BW=IB3-1                                                                  
      BMOI3=BMOI2+(BW*PAS4)                                                     
      BPLU3=BMOI3+PAS4                                                          
      KRBM3=COEFF/(BMOI3*BMOI3)                                                 
      KRBP3=COEFF/(BPLU3*BPLU3)                                                 
      DO 4200 IC3=1,2                                                           
      CW=IC3-1                                                                  
      CMOI3=CMOI2+(CW*PAS4)                                                     
      ATEST=ABS(AMOI3-AA)                                                       
      BTEST=ABS(BMOI3-BB)                                                       
      CTEST=ABS(CMOI3-CC)                                                       
      IF(ATEST.LT.0.25.AND.BTEST.LT.0.25.AND.CTEST.LT.0.25) GO TO 4200          
      CPLU3=CMOI3+PAS4                                                          
      VTESTM=AMOI3*BMOI3*CMOI3                                                  
      IF(VSUP.LT.VTESTM)GO TO 4100                                              
      VTESTP=APLU3*BPLU3*CPLU3                                                  
      IF(VINF.GT.VTESTP)GO TO 4200                                              
      KRCM3=COEFF/(CMOI3*CMOI3)                                                 
      KRCP3=COEFF/(CPLU3*CPLU3)                                                 
      CALL ORDHKL(KRAM3,KRAP3,KRBM3,KRBP3,KRCM3,KRCP3)                          
      DO 122 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 4200                                                 
  122 CONTINUE                                                                  
      NDICH(3)=NDICH(3)+1                                                       
      DO 5000 IA4=1,2                                                           
      AW=IA4-1                                                                  
      AMOI4=AMOI3+(AW*PAS8)                                                     
      APLU4=AMOI4+PAS8                                                          
      KRAM4=COEFF/(AMOI4*AMOI4)                                                 
      KRAP4=COEFF/(APLU4*APLU4)                                                 
      DO 5100 IB4=1,2                                                           
      BW=IB4-1                                                                  
      BMOI4=BMOI3+(BW*PAS8)                                                     
      BPLU4=BMOI4+PAS8                                                          
      KRBM4=COEFF/(BMOI4*BMOI4)                                                 
      KRBP4=COEFF/(BPLU4*BPLU4)                                                 
      DO 5200 IC4=1,2                                                           
      CW=IC4-1                                                                  
      CMOI4=CMOI3+(CW*PAS8)                                                     
      ATEST=ABS(AMOI4-AA)                                                       
      BTEST=ABS(BMOI4-BB)                                                       
      CTEST=ABS(CMOI4-CC)                                                       
      IF(ATEST.LT.0.25.AND.BTEST.LT.0.25.AND.CTEST.LT.0.25) GO TO 5200          
      CPLU4=CMOI4+PAS8                                                          
      VTESTM=AMOI4*BMOI4*CMOI4                                                  
      IF(VSUP.LT.VTESTM)GO TO 5100                                              
      VTESTP=APLU4*BPLU4*CPLU4                                                  
      IF(VINF.GT.VTESTP)GO TO 5200                                              
      KRCM4=COEFF/(CMOI4*CMOI4)                                                 
      KRCP4=COEFF/(CPLU4*CPLU4)                                                 
      CALL ORDHKL(KRAM4,KRAP4,KRBM4,KRBP4,KRCM4,KRCP4)                          
      DO 123 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 5200                                                 
  123 CONTINUE                                                                  
      NDICH(4)=NDICH(4)+1                                                       
      DO 6000 IA5=1,2                                                           
      AW=IA5-1                                                                  
      AMOI5=AMOI4+(AW*PAS16)                                                    
      APLU5=AMOI5+PAS16                                                         
      KRAM5=COEFF/(AMOI5*AMOI5)                                                 
      KRAP5=COEFF/(APLU5*APLU5)                                                 
      DO 6100 IB5=1,2                                                           
      BW=IB5-1                                                                  
      BMOI5=BMOI4+(BW*PAS16)                                                    
      IF(BMOI5.GT.APLU5)GO TO 6000                                              
      BPLU5=BMOI5+PAS16                                                         
      KRBM5=COEFF/(BMOI5*BMOI5)                                                 
      KRBP5=COEFF/(BPLU5*BPLU5)                                                 
      DO 6200 IC5=1,2                                                           
      CW=IC5-1                                                                  
      CMOI5=CMOI4+(CW*PAS16)                                                    
      IF(CMOI5.GT.BPLU5)GO TO 6100                                              
      ATEST=ABS(AMOI5-AA)                                                       
      BTEST=ABS(BMOI5-BB)                                                       
      CTEST=ABS(CMOI5-CC)                                                       
      IF(ATEST.LT.0.25.AND.BTEST.LT.0.25.AND.CTEST.LT.0.25) GO TO 6200          
      CPLU5=CMOI5+PAS16                                                         
      VTESTM=AMOI5*BMOI5*CMOI5                                                  
      IF(VSUP.LT.VTESTM)GO TO 6100                                              
      VTESTP=APLU5*BPLU5*CPLU5                                                  
      IF(VINF.GT.VTESTP)GO TO 6200                                              
      KRCM5=COEFF/(CMOI5*CMOI5)                                                 
      KRCP5=COEFF/(CPLU5*CPLU5)                                                 
      CALL ORDHKL(KRAM5,KRAP5,KRBM5,KRBP5,KRCM5,KRCP5)                          
      DO 124 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 6200                                                 
  124 CONTINUE                                                                  
      NDICH(5)=NDICH(5)+1                                                       
      DO 7000 IA6=1,2                                                           
      AW=IA6-1                                                                  
      AMOI6=AMOI5+(AW*PAS32)                                                    
      APLU6=AMOI6+PAS32                                                         
      KRAM6=COEFF/(AMOI6*AMOI6)                                                 
      KRAP6=COEFF/(APLU6*APLU6)                                                 
      DO 7100 IB6=1,2                                                           
      BW=IB6-1                                                                  
      BMOI6=BMOI5+(BW*PAS32)                                                    
      BPLU6=BMOI6+PAS32                                                         
      KRBM6=COEFF/(BMOI6*BMOI6)                                                 
      KRBP6=COEFF/(BPLU6*BPLU6)                                                 
      DO 7200 IC6=1,2                                                           
      CW=IC6-1                                                                  
      CMOI6=CMOI5+(CW*PAS32)                                                    
      ATEST=ABS(AMOI6-AA)                                                       
      BTEST=ABS(BMOI6-BB)                                                       
      CTEST=ABS(CMOI6-CC)                                                       
      IF(ATEST.LT.0.25.AND.BTEST.LT.0.25.AND.CTEST.LT.0.25) GO TO 7200          
      CPLU6=CMOI6+PAS32                                                         
      VTESTM=AMOI6*BMOI6*CMOI6                                                  
      IF(VSUP.LT.VTESTM)GO TO 7100                                              
      VTESTP=APLU6*BPLU6*CPLU6                                                  
      IF(VINF.GT.VTESTP)GO TO 7200                                              
      KRCM6=COEFF/(CMOI6*CMOI6)                                                 
      KRCP6=COEFF/(CPLU6*CPLU6)                                                 
      CALL ORDHKL(KRAM6,KRAP6,KRBM6,KRBP6,KRCM6,KRCP6)                          
      DO 125 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 7200                                                 
  125 CONTINUE                                                                  
      NDICH(6)=NDICH(6)+1                                                       
      DO 8000 IA7=1,2                                                           
      AW=IA7-1                                                                  
      AMOI7=AMOI6+(AW*PAS64)                                                    
      APLU7=AMOI7+PAS64                                                         
      RAM7=1./(AMOI7*AMOI7)                                                     
      RAP7=1./(APLU7*APLU7)                                                     
      KRAM7=COEFF*RAM7                                                          
      KRAP7=COEFF*RAP7                                                          
      DO 8100 IB7=1,2                                                           
      BW=IB7-1                                                                  
      BMOI7=BMOI6+(BW*PAS64)                                                    
      BPLU7=BMOI7+PAS64                                                         
      RBM7=1./(BMOI7*BMOI7)                                                     
      RBP7=1./(BPLU7*BPLU7)                                                     
      KRBM7=COEFF*RBM7                                                          
      KRBP7=COEFF*RBP7                                                          
      DO 8200 IC7=1,2                                                           
      CW=IC7-1                                                                  
      CMOI7=CMOI6+(CW*PAS64)                                                    
      ATEST=ABS(AMOI7-AA)                                                       
      BTEST=ABS(BMOI7-BB)                                                       
      CTEST=ABS(CMOI7-CC)                                                       
      IF(ATEST.LT.0.25.AND.BTEST.LT.0.25.AND.CTEST.LT.0.25)GO TO 8200           
      CPLU7=CMOI7+PAS64                                                         
      VTESTM=AMOI7*BMOI7*CMOI7                                                  
      IF(VSUP.LT.VTESTM)GO TO 8100                                              
      VTESTP=APLU7*BPLU7*CPLU7                                                  
      IF(VINF.GT.VTESTP)GO TO 8200                                              
      RCM7=1./(CMOI7*CMOI7)                                                     
      RCP7=1./(CPLU7*CPLU7)                                                     
      KRCM7=COEFF*RCM7                                                          
      KRCP7=COEFF*RCP7                                                          
      CALL ORDHKL(KRAM7,KRAP7,KRBM7,KRBP7,KRCM7,KRCP7)                          
      DO 126 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 8200                                                 
  126 CONTINUE                                                                  
      NDICH(7)=NDICH(7)+1                                                       
      AA=AMOI7+(PAS64/2.)                                                       
      BB=BMOI7+(PAS64/2.)                                                       
      CC=CMOI7+(PAS64/2.)                                                       
      BETA=90.                                                                  
      NPOSS=0                                                                   
      DO 15 I=1,N                                                               
      DO 16 J=1,10                                                              
   16 QCAL(I,J)=-1.                                                             
   15 CONTINUE                                                                  
      DO 80 M=1,MH2                                                             
      MXH=M-1                                                                   
      CARH=MXH*MXH                                                              
      DM1=CARH*RAM7                                                             
      DP1=CARH*RAP7                                                             
      DO 90 K=1,MK2                                                             
      MXK=K-1                                                                   
      CARK=MXK*MXK                                                              
      DM2=DM1+(CARK*RBM7)                                                       
      DP2=DP1+(CARK*RBP7)                                                       
      DO 1001 L=1,ML2                                                           
      MXL=L-1                                                                   
      CARL=MXL*MXL                                                              
      QQM=DP2+(CARL*RCP7)                                                       
      QQP=DM2+(CARL*RCM7)                                                       
      BORINF=QQM-EPSQ(N)                                                        
      IF(BORINF.GT.Q(N))GO TO 90                                                
      BORSUP=QQP+EPSQ(N)                                                        
      IF (BORSUP.LE.Q(N)) NPOSS=NPOSS+1                                         
      DO 110 I=1,N                                                              
      IF(QQP.LT.(Q(I)-EPSQ(I)))GO TO 110                                        
      IF(QQM.GT.(Q(I)+EPSQ(I)))GO TO 110                                        
      J=1                                                                       
  111 IF(QCAL(I,J)) 112,113,113                                                 
  113 J=J+1                                                                     
      IF(J-10)111,111,110                                                       
  112 QCAL(I,J)=(CARH/AA**2)+(CARK/BB**2)+(CARL/CC**2)                          
      XH(I,J)=MXH                                                               
      XK(I,J)=MXK                                                               
      XL(I,J)=MXL                                                               
  110 CONTINUE                                                                  
 1001 CONTINUE                                                                  
   90 CONTINUE                                                                  
   80 CONTINUE                                                                  
      DO 131 I=1,N                                                              
      IF(QCAL(I,1))8200,131,131                                                 
  131 CONTINUE                                                                  
      IND=4                                                                     
      NRIND=3                                                                   
      CALL AFFPAR(IND,NRIND,NPOSS)                                              
      WRITE(6,3)                                                                
 8200 CONTINUE                                                                  
 8100 CONTINUE                                                                  
 8000 CONTINUE                                                                  
 7200 CONTINUE                                                                  
 7100 CONTINUE                                                                  
 7000 CONTINUE                                                                  
 6200 CONTINUE                                                                  
 6100 CONTINUE                                                                  
 6000 CONTINUE                                                                  
 5200 CONTINUE                                                                  
 5100 CONTINUE                                                                  
 5000 CONTINUE                                                                  
 4200 CONTINUE                                                                  
 4100 CONTINUE                                                                  
 4000 CONTINUE                                                                  
 3200 CONTINUE
      IC2=IC2+1
      IF(IC2.LE.2) GO TO 3210                                              
 3100 CONTINUE
      IB2=IB2+1
      IF(IB2.LE.2) GO TO 3110                                           
 3000 CONTINUE
      IA2=IA2+1
      IF(IA2.LE.2) GO TO 3010                                                   
 2500 CONTINUE
      IC=IC+1
      IF(IC.LE.NC) GO TO 2510                                                   
 2000 CONTINUE
      IB=IB+1
      IF(IB.LE.NB) GO TO 2010                                                   
 1000 CONTINUE
      IA=IA+1
      IF(IA.LE.NA) GO TO 1010                                                   
      WRITE(6,10)(NDICH(I),I=1,7)                                               
   10 FORMAT(1H0,10X,42HNOMBRE D ITERATIONS DANS LES 7 DICHOTOMIES/1H ,1        
     10X,7(3X,I10))                                                             
      RETURN                                                                    
      END                                                                       
      SUBROUTINE ORDHKL(KRAM,KRAP,KRBM,KRBP,KRCM,KRCP)                          
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      DO 15 I=1,N                                                               
      KQT(I)=-1                                                                 
   15 CONTINUE                                                                  
      DO 10 M=1,MH                                                              
      MCARH=(M-1)*(M-1)                                                         
      M1=MCARH*KRAM                                                             
      MP1=MCARH*KRAP                                                            
      DO 20 K=1,MK                                                              
      MCARK=(K-1)*(K-1)                                                         
      M2=(MCARK*KRBM)+M1                                                        
      MP2=(MCARK*KRBP)+MP1                                                      
      KBORIN=MP2-KEPSQ(5)                                                       
      IF(KBORIN.GT.KQ(5))GO TO 10                                               
      DO 30 L=1,ML                                                              
      MCARL=(L-1)*(L-1)                                                         
      KQQM=(MCARL*KRCP)+MP2                                                     
      KBORIN=KQQM-KEPSQ(5)                                                      
      IF(KBORIN.GT.KQ(5))GO TO 20                                               
      KQQP=(MCARL*KRCM)+M2                                                      
      KBORSU=KQQP+KEPSQ(1)                                                      
      IF(KQ(1).GT.KBORSU)GO TO 30                                               
      DO 40 I=1,5                                                               
      KBORIN=KQQM-KEPSQ(I)                                                      
      IF(KBORIN.GT.KQ(I))GO TO 40                                               
      KBORSU=KQQP+KEPSQ(I)                                                      
      IF(KBORSU.GE.KQ(I))KQT(I)=1                                               
   40 CONTINUE                                                                  
   30 CONTINUE                                                                  
   20 CONTINUE                                                                  
   10 CONTINUE                                                                  
      DO 70 I=1,5                                                               
      IF(KQT(I).LT.0)GO TO 3000                                                 
   70 CONTINUE                                                                  
      DO 80 M=1,MH2                                                             
      MCARH=(M-1)*(M-1)                                                         
      M1=MCARH*KRAM                                                             
      MP1=MCARH*KRAP                                                            
      DO 90 K=1,MK2                                                             
      MCARK=(K-1)*(K-1)                                                         
      M2=(MCARK*KRBM)+M1                                                        
      MP2=(MCARK*KRBP)+MP1                                                      
      KBORIN=MP2-KEPSQ(N)                                                       
      IF(KBORIN.GT.KQ(N))GO TO 80                                               
      DO 100 L=1,ML2                                                            
      MCARL=(L-1)*(L-1)                                                         
      KQQP=(MCARL*KRCM)+M2                                                      
      KBORSU=KQQP+KEPSQ(6)                                                      
      IF(KBORSU.LT.KQ(6))GO TO 100                                              
      KQQM=(MCARL*KRCP)+MP2                                                     
      KBORIN=KQQM-KEPSQ(N)                                                      
      IF(KBORIN.GT.KQ(N))GO TO 90                                               
      DO 110 I=6,N                                                              
      KBORIN=KQQM-KEPSQ(I)                                                      
      IF(KBORIN.GT.KQ(I))GO TO 110                                              
      KBORSU=KQQP+KEPSQ(I)                                                      
      IF(KBORSU.GE.KQ(I))KQT(I)=1                                               
  110 CONTINUE                                                                  
  100 CONTINUE                                                                  
   90 CONTINUE                                                                  
   80 CONTINUE                                                                  
 3000 RETURN                                                                    
      END                                                                       
C     ------------------------------------------------------------------        
C     !                    M O N O C L I N I Q U E                     !        
C     ------------------------------------------------------------------        
      SUBROUTINE MONOC(AMIN,NA,NB,CMIN,NC,BEMAX,BEMIN)                          
      DIMENSION NDICH(7)                                                        
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      COMMON QCAL(50,11),XH(50,11),XK(50,11),XL(50,11),AA,BB,CC,BETA            
    1 FORMAT(1H0,5X,2HA=,F8.5,5X,2HB=,F8.5,5X,2HC=,F8.5,5X,5HBETA=,F8.3,        
     15X,33HVOLUME DE LA MAILLE ELEMENTAIRE =,F10.3)                            
    2 FORMAT(1H ,10X,F8.5,4X,3(10X,F3.0),9X,F8.5,11X,F8.5,9X,F8.5,10X,F8        
     1.5)                                                                       
    3 FORMAT(1H0,50X,9(1H-)/)                                                   
    4 FORMAT(1H ,20X,23HDOMAINE ANGLE EXAMINE :,5X,13HBETA MINIMUM=,F10.        
     12,7H DEGRES,5X,13HBETA MAXIMUM=,F10.2,7H DEGRES)                
    5 FORMAT(1H0,5X,2HA=,F8.5,5X,2HB=,F8.5,5X,2HC=,F8.5,5X,5HBETA=,F8.3,        
     15X,33HVOLUME DE LA MAILLE ELEMENTAIRE =,F10.3,4X,2HZ=,I2)                 
    6 FORMAT(1H0,7X,14HD EXPERIMENTAL,12X,1HH,12X,1HK,12X,1HL,10X,9HD CA        
     1LCULE,12X,5HDE-DC,12X,7HEPSILON/1H0)                                      
      DO 20 I=1,7                                                               
   20 NDICH(I)=0                                                                
      AA=100.                                                                   
      BB=100.                                                                   
      CC=100.                                                                   
      BEBE=3.                                                                   
      BMIN=CMIN                                                                 
      PASBE=PIRAD*2.5                                                           
      NBETA=((BEMAX-BEMIN)/PASBE)+1.                                            
      PAB2=PASBE/2.                                                             
      PAB4=PASBE/4.                                                             
      PAB8=PASBE/8.                                                             
      PAB16=PASBE/16.                                                           
      PAB32=PASBE/32.                                                           
      PAB64=PASBE/64.                                                           
      DO 1000 IBET1=1,NBETA                                                     
      BETW=IBET1-1                                                              
      BEINF1=BEMIN+(BETW*PASBE)                                                 
      BESUP1=BEINF1+PASBE                                                       
      BBETM=BEINF1*PIDEG                                                        
      BBETP=BESUP1*PIDEG                                                        
      WRITE(6,4)BBETM,BBETP 
      WRITE(85,4)BBETM,BBETP 
      SM1=1./(SIN(BEINF1))**2                                                   
      SP1=1./(SIN(BESUP1))**2                                                   
      CSM1=COS(BEINF1)*SM1                                                      
      CSP1=COS(BESUP1)*SP1                                                      
      DO 1100 IA=1,NA                                                           
      AW=IA-1                                                                   
      AMOI1=AMIN+(AW*PAS)                                                       
      APLU1=AMOI1+PAS                                                           
      APMO=APLU1*SIN(BEINF1)                                                    
      DO 21 I=1,10                                                              
      PP=I                                                                      
      DA=APMO/PP                                                                
      IF(D(5)-DA)21,22,22                                                       
   21 CONTINUE                                                                  
   22 MH=I                                                                      
      DO 55 J=I,20                                                              
      PP=J                                                                      
      DA=APMO/PP                                                                
      IF(D(N)-DA)55,51,51                                                       
   55 CONTINUE                                                                  
   51 MH2=J                                                                     
      KRAP1=COEFF*SP1/(AMOI1*AMOI1)                                             
      KRAM1=COEFF*SM1/(APLU1*APLU1)                                             
      DO 1200 IB=1,NB                                                           
      BW=IB-1                                                                   
      BMOI1=BMIN+(BW*PAS)                                                       
      BPLU1=BMOI1+PAS                                                           
      IF((APLU1.LT.D(1)).AND.(BPLU1.LT.D(1)))GO TO 1200                         
      DO 23 I=1,10                                                              
      PP=I                                                                      
      DB=BPLU1/PP                                                               
      IF(D(5)-DB)23,24,24                                                       
   23 CONTINUE                                                                  
   24 MK=I                                                                      
      DO 85 J=I,20                                                              
      PP=J                                                                      
      DB=BPLU1/PP                                                               
      IF(D(N)-DB)85,61,61                                                       
   85 CONTINUE                                                                  
   61 MK2=J                                                                     
      KRBM1=COEFF/(BPLU1*BPLU1)                                                 
      KRBP1=COEFF/(BMOI1*BMOI1)                                                 
      DO 1300 IC=1,NC                                                           
      CW=IC-1                                                                   
      CMOI1=CMIN+(CW*PAS)                                                       
      CPLU1=CMOI1+PAS                                                           
      IF(CMOI1.GT.APLU1)GO TO 1200                                              
      VTESTM=AMOI1*BMOI1*CMOI1*SIN(BESUP1)                                      
      IF(V.LT.VTESTM)GO TO 1200                                                 
      IF(VSUP.LT.VTESTM)GO TO 1200                                              
      VTESTP=APLU1*BPLU1*CPLU1*SIN(BEINF1)                                      
      IF(VINF.GT.VTESTP)GO TO 1300                                              
      CPMO=CPLU1*SIN(BEINF1)                                                    
      DO 25 I=1,10                                                              
      PP=I                                                                      
      DC=CPMO/PP                                                                
      IF(D(5)-DC)25,26,26                                                       
   25 CONTINUE                                                                  
   26 ML=I                                                                      
      DO 75 J=I,20                                                              
      PP=J                                                                      
      DC=CPMO/PP                                                                
      IF(D(N)-DC)75,71,71                                                       
   75 CONTINUE                                                                  
   71 ML2=J                                                                     
      KRCM=COEFF*SM1/(CPLU1*CPLU1)                                              
      KRCP=COEFF*SP1/(CMOI1*CMOI1)                                              
      KRACM=2.*COEFF*CSM1/(APLU1*CPLU1)                                         
      KRACP=2.*COEFF*CSP1/(AMOI1*CMOI1)                                         
      CALL MODHKL(KRAM1,KRAP1,KRCM,KRCP,KRACM,KRACP,KRBM1,KRBP1)                
      DO 120 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 1300                                                 
  120 CONTINUE                                                                  
      NDICH(1)=NDICH(1)+1                                                       
      DO 2000 IBET2=1,2                                                         
      BETW=IBET2-1                                                              
      BEINF2=BEINF1+(BETW*PAB2)                                                 
      BESUP2=BEINF2+PAB2                                                        
      SM2=1./(SIN(BEINF2))**2                                                   
      SP2=1./(SIN(BESUP2))**2                                                   
      CSM2=COS(BEINF2)*SM2                                                      
      CSP2=COS(BESUP2)*SP2                                                      
      DO 2100 IA2=1,2                                                           
      AW=IA2-1                                                                  
      AMOI2=AMOI1+(AW*PAS2)                                                     
      APLU2=AMOI2+PAS2                                                          
      KRAM2=COEFF*SM2/(APLU2*APLU2)                                             
      KRAP2=COEFF*SP2/(AMOI2*AMOI2)                                             
      DO 2200 IB2=1,2                                                           
      BW=IB2-1                                                                  
      BMOI2=BMOI1+(BW*PAS2)                                                     
      BPLU2=BMOI2+PAS2                                                          
      KRBM2=COEFF/(BPLU2*BPLU2)                                                 
      KRBP2=COEFF/(BMOI2*BMOI2)                                                 
      DO 2300 IC2=1,2                                                           
      CW=IC2-1                                                                  
      CMOI2=CMOI1+(CW*PAS2)                                                     
      ATEST=ABS(AMOI2-AA)                                                       
      BTEST=ABS(BMOI2-BB)                                                       
      CTEST=ABS(CMOI2-CC)                                                       
      BETEST=ABS(BEINF2-BEBE)                                                   
      IF(ATEST.LT.0.1.AND.BTEST.LT.0.1.AND.CTEST.LT.0.1.AND.BETEST.LT.0.        
     1009)GO TO 2300                                                            
      VTESTM=AMOI2*BMOI2*CMOI2*SIN(BESUP2)                                      
      IF(VTESTM.GT.VSUP)GO TO 2200                                              
      CPLU2=CMOI2+PAS2                                                          
      VTESTP=APLU2*BPLU2*CPLU2*SIN(BEINF2)                                      
      IF(VINF.GT.VTESTP)GO TO 2300                                              
      KRCM=COEFF*SM2/(CPLU2*CPLU2)                                              
      KRCP=COEFF*SP2/(CMOI2*CMOI2)                                              
      KRACM=2.*COEFF*CSM2/(APLU2*CPLU2)                                         
      KRACP=2.*COEFF*CSP2/(AMOI2*CMOI2)                                         
      CALL MODHKL(KRAM2,KRAP2,KRCM,KRCP,KRACM,KRACP,KRBM2,KRBP2)                
      DO 122 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 2300                                                 
  122 CONTINUE                                                                  
      NDICH(2)=NDICH(2)+1                                                       
      DO 3000 IBET3=1,2                                                         
      BETW=IBET3-1                                                              
      BEINF3=BEINF2+(BETW*PAB4)                                                 
      BESUP3=BEINF3+PAB4                                                        
      SM3=1./(SIN(BEINF3))**2                                                   
      SP3=1./(SIN(BESUP3))**2                                                   
      CSM3=COS(BEINF3)*SM3                                                      
      CSP3=COS(BESUP3)*SP3                                                      
      DO 3100 IA3=1,2                                                           
      AW=IA3-1                                                                  
      AMOI3=AMOI2+(AW*PAS4)                                                     
      APLU3=AMOI3+PAS4                                                          
      KRAM3=COEFF*SM3/(APLU3*APLU3)                                             
      KRAP3=COEFF*SP3/(AMOI3*AMOI3)                                             
      DO 3200 IB3=1,2                                                           
      BW=IB3-1                                                                  
      BMOI3=BMOI2+(BW*PAS4)                                                     
      BPLU3=BMOI3+PAS4                                                          
      KRBM3=COEFF/(BPLU3*BPLU3)                                                 
      KRBP3=COEFF/(BMOI3*BMOI3)                                                 
      DO 3300 IC3=1,2                                                           
      CW=IC3-1                                                                  
      CMOI3=CMOI2+(CW*PAS4)                                                     
      ATEST=ABS(AMOI3-AA)                                                       
      BTEST=ABS(BMOI3-BB)                                                       
      CTEST=ABS(CMOI3-CC)                                                       
      BETEST=ABS(BEINF3-BEBE)                                                   
      IF(ATEST.LT.0.1.AND.BTEST.LT.0.1.AND.CTEST.LT.0.1.AND.BETEST.LT.0.        
     1009)GO TO 3300                                                            
      VTESTM=AMOI3*BMOI3*CMOI3*SIN(BESUP3)                                      
      IF(VTESTM.GT.VSUP)GO TO 3200                                              
      CPLU3=CMOI3+PAS4                                                          
      VTESTP=APLU3*BPLU3*CPLU3*SIN(BEINF3)                                      
      IF(VINF.GT.VTESTP)GO TO 3300                                              
      KRCM=COEFF*SM3/(CPLU3*CPLU3)                                              
      KRCP=COEFF*SP3/(CMOI3*CMOI3)                                              
      KRACM=2.*COEFF*CSM3/(APLU3*CPLU3)                                         
      KRACP=2.*COEFF*CSP3/(AMOI3*CMOI3)                                         
      CALL MODHKL(KRAM3,KRAP3,KRCM,KRCP,KRACM,KRACP,KRBM3,KRBP3)                
      DO 124 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 3300                                                 
  124 CONTINUE                                                                  
      NDICH(3)=NDICH(3)+1                                                       
      DO 4000 IBET4=1,2                                                         
      BETW=IBET4-1                                                              
      BEINF4=BEINF3+(BETW*PAB8)                                                 
      BESUP4=BEINF4+PAB8                                                        
      SM4=1./(SIN(BEINF4))**2                                                   
      SP4=1./(SIN(BESUP4))**2                                                   
      CSM4=COS(BEINF4)*SM4                                                      
      CSP4=COS(BESUP4)*SP4                                                      
      DO 4100 IA4=1,2                                                           
      AW=IA4-1                                                                  
      AMOI4=AMOI3+(AW*PAS8)                                                     
      APLU4=AMOI4+PAS8                                                          
      KRAM4=COEFF*SM4/(APLU4*APLU4)                                             
      KRAP4=COEFF*SP4/(AMOI4*AMOI4)                                             
      DO 4200 IB4=1,2                                                           
      BW=IB4-1                                                                  
      BMOI4=BMOI3+(BW*PAS8)                                                     
      BPLU4=BMOI4+PAS8                                                          
      KRBM4=COEFF/(BPLU4*BPLU4)                                                 
      KRBP4=COEFF/(BMOI4*BMOI4)                                                 
      DO 4300 IC4=1,2                                                           
      CW=IC4-1                                                                  
      CMOI4=CMOI3+(CW*PAS8)                                                     
      ATEST=ABS(AMOI4-AA)                                                       
      BTEST=ABS(BMOI4-BB)                                                       
      CTEST=ABS(CMOI4-CC)                                                       
      BETEST=ABS(BEINF4-BEBE)                                                   
      IF(ATEST.LT.0.1.AND.BTEST.LT.0.1.AND.CTEST.LT.0.1.AND.BETEST.LT.0.        
     1009)GO TO 4300                                                            
      VTESTM=AMOI4*BMOI4*CMOI4*SIN(BESUP4)                                      
      IF(VTESTM.GT.VSUP)GO TO 4200                                              
      CPLU4=CMOI4+PAS8                                                          
      VTESTP=APLU4*BPLU4*CPLU4*SIN(BEINF4)                                      
      IF(VINF.GT.VTESTP)GO TO 4300                                              
      KRCM=COEFF*SM4/(CPLU4*CPLU4)                                              
      KRCP=COEFF*SP4/(CMOI4*CMOI4)                                              
      KRACM=2.*COEFF*CSM4/(APLU4*CPLU4)                                         
      KRACP=2.*COEFF*CSP4/(AMOI4*CMOI4)                                         
      CALL MODHKL(KRAM4,KRAP4,KRCM,KRCP,KRACM,KRACP,KRBM4,KRBP4)                
      DO 125 I=6,N                                                              
      IF(KQT(I).LT.0)GO TO 4300                                                 
  125 CONTINUE                                                                  
      NDICH(4)=NDICH(4)+1                                                       
      DO 5000 IBET5=1,2                                                         
      BETW=IBET5-1                                                              
      BEINF5=BEINF4+(BETW*PAB16)                                                
      BESUP5=BEINF5+PAB16                                                       
      SM5=1./(SIN(BEINF5))**2                                                   
      SP5=1./(SIN(BESUP5))**2                                                   
      CSM5=COS(BEINF5)*SM5                                                      
      CSP5=COS(BESUP5)*SP5                                                      
      DO 5100 IA5=1,2                                                           
      AW=IA5-1                                                                  
      AMOI5=AMOI4+(AW*PAS16)                                                    
      APLU5=AMOI5+PAS16                                                         
      KRAM5=COEFF*SM5/(APLU5*APLU5)                                             
      KRAP5=COEFF*SP5/(AMOI5*AMOI5)                                             
      DO 5200 IB5=1,2                                                           
      BW=IB5-1                                                                  
      BMOI5=BMOI4+(BW*PAS16)                                                    
      BPLU5=BMOI5+PAS16                                                         
      KRBM5=COEFF/(BPLU5*BPLU5)                                                 
      KRBP5=COEFF/(BMOI5*BMOI5)                                                 
      DO 5300 IC5=1,2                                                           
      CW=IC5-1                                                                  
      CMOI5=CMOI4+(CW*PAS16)                                                    
      ATEST=ABS(AMOI5-AA)                                                       
      BTEST=ABS(BMOI5-BB)                                                       
      CTEST=ABS(CMOI5-CC)                                                       
      BETEST=ABS(BEINF5-BEBE)                                                   
      IF(ATEST.LT.0.1.AND.BTEST.LT.0.1.AND.CTEST.LT.0.1.AND.BETEST.LT.0.        
     1009)GO TO 5300                                                            
      VTESTM=AMOI5*BMOI5*CMOI5*SIN(BESUP5)                                      
      IF(VTESTM.GT.VSUP)GO TO 5200                                              
      CPLU5=CMOI5+PAS16                                                         
      VTESTP=APLU5*BPLU5*CPLU5*SIN(BEINF5)                                      
      IF(VINF.GT.VTESTP)GO TO 5300                                              
      KRCM=COEFF*SM5/(CPLU5*CPLU5)                                              
      KRCP=COEFF*SP5/(CMOI5*CMOI5)                                              
      KRACM=2.*COEFF*CSM5/(APLU5*CPLU5)                                         
      KRACP=2.*COEFF*CSP5/(AMOI5*CMOI5)                                         
      CALL MODHKL(KRAM5,KRAP5,KRCM,KRCP,KRACM,KRACP,KRBM5,KRBP5)                
CCCC      DO 126 I=6,N
      I=6                                                              
126      IF(KQT(I).LT.0)GO TO 5300                                       
CCCC  126 CONTINUE
      I=I+1
      IF(I.LE.N) GO TO 126                                                    
      NDICH(5)=NDICH(5)+1                                                       
CCCC      DO 6000 IBET6 =1,2
      IBET6=1                                                        
6001      BETW=IBET6-1                                                  
      BEINF6=BEINF5+(BETW*PAB32)                                                
      BESUP6=BEINF6+PAB32                                                       
      SM6=1./(SIN(BEINF6))**2                                                   
      SP6=1./(SIN(BESUP6))**2                                                   
      CSM6=COS(BEINF6)*SM6                                                      
      CSP6=COS(BESUP6)*SP6                                                      
CCCC      DO 6100 IA6=1,2
      IA6=1                                                           
6101      AW=IA6-1                                                       
      AMOI6=AMOI5+(AW*PAS32)                                                    
      APLU6=AMOI6+PAS32                                                         
      KRAM6=COEFF*SM6/(APLU6*APLU6)                                             
      KRAP6=COEFF*SP6/(AMOI6*AMOI6)                                             
CCCC      DO 6200 IB6=1,2
      IB6=1                                                           
6201      BW=IB6-1                                                        
      BMOI6=BMOI5+(BW*PAS32)                                                    
      BPLU6=BMOI6+PAS32                                                         
      KRBM6=COEFF/(BPLU6*BPLU6)                                                 
      KRBP6=COEFF/(BMOI6*BMOI6)                                                 
CCCC      DO 6300 IC6=1,2
      IC6=1                                                           
6301      CW=IC6-1                                                           
      CMOI6=CMOI5+(CW*PAS32)                                                    
      ATEST=ABS(AMOI6-AA)                                                       
      BTEST=ABS(BMOI6-BB)                                                       
      CTEST=ABS(CMOI6-CC)                                                       
      BETEST=ABS(BEINF6-BEBE)                                                   
      IF(ATEST.LT.0.1.AND.BTEST.LT.0.1.AND.CTEST.LT.0.1.AND.BETEST.LT.0.        
     1009)GO TO 6300                                                            
      VTESTM=AMOI6*BMOI6*CMOI6*SIN(BESUP6)                                      
      IF(VTESTM.GT.VSUP)GO TO 6200                                              
      CPLU6=CMOI6+PAS32                                                         
      VTESTP=APLU6*BPLU6*CPLU6*SIN(BEINF6)                                      
      IF(VINF.GT.VTESTP)GO TO 6300                                              
      KRCM=COEFF*SM6/(CPLU6*CPLU6)                                              
      KRCP=COEFF*SP6/(CMOI6*CMOI6)                                              
      KRACM=2.*COEFF*CSM6/(APLU6*CPLU6)                                         
      KRACP=2.*COEFF*CSP6/(AMOI6*CMOI6)                                         
      CALL MODHKL(KRAM6,KRAP6,KRCM,KRCP,KRACM,KRACP,KRBM6,KRBP6)                
CCCC      DO 127 I=6,N
      I=6                                                              
127      IF(KQT(I).LT.0)GO TO 6300                                 
CCCC  127 CONTINUE
      I=I+1
      IF(I.LE.N) GO TO 127                                                   
      NDICH(6)=NDICH(6)+1                                                       
CCCC      DO 7000 IBET7=1,2
      IBET7=1                                                         
7001      BETW=IBET7-1                                                     
      BEINF7=BEINF6+(BETW*PAB64)                                                
      BESUP7=BEINF7+PAB64                                                       
      SM7=1./(SIN(BEINF7))**2                                                   
      SP7=1./(SIN(BESUP7))**2                                                   
      CSM7=COS(BEINF7)*SM7                                                      
      CSP7=COS(BESUP7)*SP7                                                      
CCCC      DO 7100 IA7=1,2
      IA7=1                                                           
7101      AW=IA7-1                                                        
      AMOI7=AMOI6+(AW*PAS64)                                                    
      APLU7=AMOI7+PAS64                                                         
      RAM7=SM7/APLU7**2                                                         
      RAP7=SP7/AMOI7**2                                                         
      KRAM7=COEFF*SM7/(APLU7*APLU7)                                             
      KRAP7=COEFF*SP7/(AMOI7*AMOI7)                                             
CCCC      DO 7200 IB7=1,2
      IB7=1                                                           
7201      BW=IB7-1                                                      
      BMOI7=BMOI6+(BW*PAS64)                                                    
      BPLU7=BMOI7+PAS64                                                         
      RBM7=1./BPLU7**2                                                          
      RBP7=1./BMOI7**2                                                          
      KRBM7=COEFF/(BPLU7*BPLU7)                                                 
      KRBP7=COEFF/(BMOI7*BMOI7)                                                 
CCCC      DO 7300 IC7=1,2
      IC7=1                                                           
7301      CW=IC7-1                                                        
      CMOI7=CMOI6+(CW*PAS64)                                                    
      ATEST=ABS(AMOI7-AA)                                                       
      BTEST=ABS(BMOI7-BB)                                                       
      CTEST=ABS(CMOI7-CC)                                                       
      BETEST=ABS(BEINF7-BEBE)                                                   
      IF(ATEST.LT.0.1.AND.BTEST.LT.0.1.AND.CTEST.LT.0.1.AND.BETEST.LT.0.        
     1009)GO TO 7300                                                            
      VTESTM=AMOI7*BMOI7*CMOI7*SIN(BESUP7)                                      
      IF(VTESTM.GT.VSUP)GO TO 7200                                              
      CPLU7=CMOI7+PAS64                                                         
      VTESTP=APLU7*BPLU7*CPLU7*SIN(BEINF7)                                      
      IF(VINF.GT.VTESTP)GO TO 7300                                              
      RCM=SM7/CPLU7**2                                                          
      RCP=SP7/CMOI7**2                                                          
      RACM=2.*CSM7/(APLU7*CPLU7)                                                
      RACP=2.*CSP7/(AMOI7*CMOI7)                                                
      KRCM=COEFF*SM7/(CPLU7*CPLU7)                                              
      KRCP=COEFF*SP7/(CMOI7*CMOI7)                                              
      KRACM=2.*COEFF*CSM7/(APLU7*CPLU7)                                         
      KRACP=2.*COEFF*CSP7/(AMOI7*CMOI7)                                         
      CALL MODHKL(KRAM7,KRAP7,KRCM,KRCP,KRACM,KRACP,KRBM7,KRBP7)                
CCCC      DO 128 I=6,N
      I=6                                                              
128      IF(KQT(I).LT.0)GO TO 7300                                       
CCCC  128 CONTINUE
      I=I+1
      IF(I.LE.N) GO TO 128                                                    
      NDICH(7)=NDICH(7)+1                                                       
      AA=AMOI7+(PAS64/2.)                                                       
      BB=BMOI7+(PAS64/2.)                                                       
      CC=CMOI7+(PAS64/2.)                                                       
      BEBE=BEINF7+(PAB64/2.)                                                    
      BETA=PIDEG*BEBE                                                           
      NPOSS=0                                                                   
CCCC      DO 15 I=1,N                                                         
CCCC      DO 16 J=1,11                                                        
CCCC   16 QCAL(I,J)=-1.                                                       
CCCC   15 CONTINUE 
      I=1
16    J=1
15    QCAL(I,J)=-1
      J=J+1
      IF(J.LT.12) GO TO 15
      I=I+1
      IF(I.LE.N)  GO TO 16                  
CCCC      DO 80 M=1,MH2
      M=1                                                             
8001      ML1=ML2*MOD(1,M)                                                   
      NML=ML1+ML2+1                                                             
      MXH=M-1                                                                   
      CARH=MXH*MXH                                                              
      DM1=CARH*RAM7                                                             
      DP1=CARH*RAP7                                                             
CCCC      DO 90 K=1,MK2
      K=1       
9001      MXK=K-1                                                            
      CARK=MXK*MXK                                                              
      DM2=(CARK*RBM7)+DM1                                                       
      DP2=(CARK*RBP7)+DP1                                                       
CCCC      DO100 L=1,NML
      L=1                                                             
10001      MXL=-ML1+L-1                                                   
      CARL=MXL*MXL                                                              
      PROHL=MXH*MXL                                                             
      DM3=(CARL*RCM)+DM2                                                        
      DP3=(CARL*RCP)+DP2                                                        
      IF(MXL)101,102,102                                                        
  101 QQM=DM3-(PROHL*RACP)                                                      
      IF((QQM-EPSQ(N)).GT.Q(N))GO TO 100                                        
      QQP=DP3-(PROHL*RACM)                                                      
      GO TO 103                                                                 
  102 QQM=DM3-(PROHL*RACM)                                                      
      IF((QQM-EPSQ(N)).GT.Q(N))GO TO 90                                         
      QQP=DP3-(PROHL*RACP)                                                      
  103 IF((QQP+EPSQ(1)).LT.Q(1))GO TO 100                                        
      NPOSS=NPOSS+1                                                             
CCCC      DO 119 I=1,N 
      I=1                                                             
1191      IF((QQM-EPSQ(I)).GT.Q(I))GO TO 119                                 
      IF((QQP+EPSQ(I)).LT.Q(I))GO TO 119                                        
      J=1                                                                       
  111 IF(QCAL(I,J))112,113,113                                                  
  113 J=J+1                                                                     
      IF(J-10)111,111,119                                                       
  112 QCAL(I,J)=(CARH/(AA**2*SIN(BEBE)**2))+(CARK/BB**2)+(CARL/(CC**2*SI        
     1N(BEBE)**2))-(2.*PROHL*COS(BEBE)/(AA*CC*SIN(BEBE)**2))                    
      XH(I,J)=MXH                                                               
      XK(I,J)=MXK                                                               
      XL(I,J)=MXL                                                               
  119 CONTINUE
      I=I+1
      IF(I.LE.N) GO TO 1191
  100 CONTINUE
      L=L+1
      IF(L.LE.NML) GO TO 10001
   90 CONTINUE
      K=K+1
      IF(K.LE.MK2) GO TO 9001
   80 CONTINUE
      M=M+1
      IF(M.LE.MH2) GO TO 8001
      IND=5                                                                     
      NRIND=4                                                                   
      CALL AFFPAR(IND,NRIND,NPOSS)                                              
      WRITE(6,3)                                                                
 7300 CONTINUE
      IC7=IC7+1
      IF(IC7.LE.2) GO TO 7301
 7200 CONTINUE
      IB7=IB7+1
      IF(IB7.LE.2) GO TO 7201                                                   
 7100 CONTINUE
      IA7=IA7+1
      IF(IA7.LE.2) GO TO 7101                                                   
 7000 CONTINUE
      IBET7=IBET7+1
      IF(IBET7.LE.2) GO TO 7001                                          
 6300 CONTINUE
      IC6=IC6+1
      IF(IC6.LE.2) GO TO 6301                                               
 6200 CONTINUE
      IB6=IB6+1
      IF(IB6.LE.2) GO TO 6201
 6100 CONTINUE
      IA6=IA6+1
      IF(IA6.LE.2) GO TO 6101                                                   
 6000 CONTINUE
      IBET6=IBET6+1
      IF(IBET6.LE.2) GO TO 6001                                           
 5300 CONTINUE                                                                  
 5200 CONTINUE                                                                  
 5100 CONTINUE                                                                  
 5000 CONTINUE                                                                  
 4300 CONTINUE                                                                  
 4200 CONTINUE                                                                  
 4100 CONTINUE                                                                  
 4000 CONTINUE                                                                  
 3300 CONTINUE                                                                  
 3200 CONTINUE                                                                  
 3100 CONTINUE                                                                  
 3000 CONTINUE                                                                  
 2300 CONTINUE                                                                  
 2200 CONTINUE                                                                  
 2100 CONTINUE                                                                  
 2000 CONTINUE                                                                  
 1300 CONTINUE                                                                  
 1200 CONTINUE                                                                  
 1100 CONTINUE                                                                  
 1000 CONTINUE                                                                  
      WRITE(6,10)(NDICH(I),I=1,7)                                               
   10 FORMAT(1H0,10X,42HNOMBRE D ITERATIONS DANS LES 7 DICHOTOMIES/1H ,1        
     10X,7(3X,I10))                                                             
      RETURN                                                                    
      END                                                                       
      SUBROUTINE MODHKL(KRAM,KRAP,KRCM,KRCP,KRACM,KRACP,KRBM,KRBP)              
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      DO 15 I=1,N                                                               
      KQT(I)=-1                                                                 
   15 CONTINUE                                                                  
      NML=2*ML+1                                                                
      DO 10 M=1,MH                                                              
      MXH=M-1                                                                   
      MCARH=MXH*MXH                                                             
      MM1=MCARH*KRAM                                                            
      MP1=MCARH*KRAP                                                            
      DO 20 K=1,MK                                                              
      MXK=K-1                                                                   
      MCARK=MXK*MXK                                                             
      MM2=(MCARK*KRBM)+MM1                                                      
      MP2=(MCARK*KRBP)+MP1                                                      
      DO 30 L=1,NML                                                             
      MXL=-ML+L-1                                                               
      MCARL=MXL*MXL                                                             
      MPROHL=MXH*MXL                                                            
      MM3=(MCARL*KRCM)+MM2                                                      
      MP3=(MCARL*KRCP)+MP2                                                      
      IF(MXL)50,51,51                                                           
   50 KQQM=MM3-(MPROHL*KRACP)                                                   
      IF((KQQM-KEPSQ(5)).GT.KQ(5))GO TO 30                                      
      KQQP=MP3-(MPROHL*KRACM)                                                   
      GO TO 52                                                                  
   51 KQQM=MM3-(MPROHL*KRACM)                                                   
      IF((KQQM-KEPSQ(5)).GT.KQ(5))GO TO 20                                      
      KQQP=MP3-(MPROHL*KRACP)                                                   
   52 IF((KQQP+KEPSQ(1)).LT.KQ(1))GO TO 30                                      
      DO 40 I=1,5                                                               
      KBORIN=KQQM-KEPSQ(I)                                                      
      IF(KBORIN.GT.KQ(I))GO TO 40                                               
      KBORSU=KQQP+KEPSQ(I)                                                      
      IF(KBORSU.GE.KQ(I))KQT(I)=1                                               
   40 CONTINUE                                                                  
   30 CONTINUE                                                                  
   20 CONTINUE                                                                  
   10 CONTINUE                                                                  
      DO 70 I=1,5                                                               
      IF(KQT(I).LT.0)GO TO 3000                                                 
   70 CONTINUE                                                                  
      NML2=2*ML2+1                                                              
      DO 80 M=1,MH2                                                             
      MXH=M-1                                                                   
      MCARH=MXH*MXH                                                             
      MM1=MCARH*KRAM                                                            
      MP1=MCARH*KRAP                                                            
      DO 90 K=1,MK2                                                             
      MXK=K-1                                                                   
      MCARK=MXK*MXK                                                             
      MM2=(MCARK*KRBM)+MM1                                                      
      MP2=(MCARK*KRBP)+MP1                                                      
      DO 100 L=1,NML2                                                           
      MXL=-ML2-1+L                                                              
      MCARL=MXL*MXL                                                             
      MPROHL=MXH*MXL                                                            
      MM3=(MCARL*KRCM)+MM2                                                      
      MP3=(MCARL*KRCP)+MP2                                                      
      IF(MXL)101,102,102                                                        
  101 KQQM=MM3-(MPROHL*KRACP)                                                   
      IF((KQQM-KEPSQ(N)).GT.KQ(N))GO TO 100                                     
      KQQP=MP3-(MPROHL*KRACM)                                                   
      GO TO 103                                                                 
  102 KQQM=MM3-(MPROHL*KRACM)                                                   
      IF((KQQM-KEPSQ(N)).GT.KQ(N))GO TO 90                                      
      KQQP=MP3-(MPROHL*KRACP)                                                   
  103 IF((KQQP+KEPSQ(6)).LT.KQ(6))GO TO 100                                     
      DO 110 I=6,N                                                              
      KBORIN=KQQM-KEPSQ(I)                                                      
      IF(KBORIN.GT.KQ(I))GO TO 110                                              
      KBORSU=KQQP+KEPSQ(I)                                                      
      IF(KBORSU.GE.KQ(I))KQT(I)=1                                               
  110 CONTINUE                                                                  
  100 CONTINUE                                                                  
   90 CONTINUE                                                                  
   80 CONTINUE                                                                  
 3000 RETURN                                                                    
      END                                                                       
C     ------------------------------------------------------------------        
C     !       A F F I N E M E N T   D E S   P A R A M E T R E S        !        
C     ------------------------------------------------------------------        
      SUBROUTINE AFFPAR(IND,NRIND,NPOSS)                                        
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      COMMON QCAL(50,11),XH(50,11),XK(50,11),XL(50,11),AA,BB,CC,BETA            
      COMMON REC(4),PAR(4)                                                      
      DIMENSION JH(50),JK(50),JL(50),POIDS(50),SSQ(50),SSC(50),DIR(4)           
      DIMENSION A(5,5),B(4,4),E(5),ECART(4),ERROR(4),TEM(4)
      DATA NUMANS/0/
      DO 710 I=1,3                                                              
      ERROR(I)=0.0001                                                           
  710 CONTINUE                                                                  
      ERROR(4)=0.001                                                            
      BIDON1=1.001                                                              
      BIDON2=1000./1001.                                                        
      PAR(1)=AA                                                                 
      PAR(2)=BB                                                                 
      PAR(3)=CC                                                                 
      PAR(4)=BETA                                                               
C     S E L E C T I O N   D E S   R A I E S   U N I V O Q U E S                 
      NRAIES=0                                                                  
      DO 10 I=1,N                                                               
      IF (QCAL(I,2).GE.0) GOTO10                                                
      NRAIES=NRAIES+1                                                           
      JH(NRAIES)=XH(I,1)                                                        
      JK(NRAIES)=XK(I,1)                                                        
      JL(NRAIES)=XL(I,1)                                                        
      POIDS(NRAIES)=2./EPSQ(I)                                                  
      SSQ(NRAIES)=Q(I)                                                          
   10 CONTINUE                                                                  
      RINDEP=NRAIES-IND+1                                                       
      IF (RINDEP.LE.0.0) GOTO 565                                               
C     CAMBIO SEGUN HEXAGONAL                                                    
      IF (IND.NE.3) GOTO30                                                      
      PAR(2)=CC                                                                 
      PAR(3)=AA                                                                 
      PAR(4)=120.                                                               
      DO 20 I=1,NRAIES                                                          
      INUTIL=JL(I)                                                              
      JL(I)=JK(I)                                                               
      JK(I)=INUTIL                                                              
   20 CONTINUE                                                                  
C     TEST N. 1 : NOMBRE DE RAIES                                               
   30 IF (NRAIES.LT.NRIND) GOTO 565                                             
      LLAVE=1                                                                   
C      A F F I N E M E N T                                                      
      CALL PASAJE                                                               
      DO 210 I=1,4                                                              
      PAR(I)=REC(I)                                                             
  210 CONTINUE                                                                  
      CALL PASAJE                                                               
      DO 215 I=1,4                                                              
      DIR(I)=REC(I)                                                             
  215 CONTINUE                                                                  
      ITERAT=0                                                                  
  600 ITERAT=ITERAT+1                                                           
      DO 220 J=1,5                                                              
      DO 225 K=1,5                                                              
      A(J,K)=0                                                                  
  225 CONTINUE                                                                  
  220 CONTINUE                                                                  
      XANGLE=PIRAD*PAR(4)                                                       
      COSX=COS(XANGLE)                                                          
      SINX=SIN(XANGLE)                                                          
      DO 230 I=1,NRAIES                                                         
      PARH=JH(I)*PAR(1)                                                         
      PARK=JK(I)*PAR(2)                                                         
      PARL=JL(I)*PAR(3)                                                         
      E(1)=2.*JH(I)*(PARH+PARL*COSX)*POIDS(I)                                   
      E(2)=2.*JK(I)*PARK*POIDS(I)                                               
      E(3)=2.*JL(I)*(PARL+PARH*COSX)*POIDS(I)                                   
      E(4)=0.                                                                   
      GOTO(235,240,245,255,250),IND                                             
  235 E(1)=E(1)+E(2)+E(3)                                                       
      E(2)=0.                                                                   
      E(3)=0.                                                                   
      GOTO 255                                                                  
  240 E(1)=E(1)+E(2)                                                            
      E(2)=0.                                                                   
      GOTO 255                                                                  
  245 E(1)=E(1)+E(3)                                                            
      E(3)=0.                                                                   
      GOTO 255                                                                  
  250 E(4)=-2.*PIRAD*PARH*PARL*SINX*POIDS(I)                                    
  255 SSC(I)=(PARH*PARH)+(PARK*PARK)+(PARL*PARL)+(2.*PARH*PARL*COSX)            
      E(5)=(SSQ(I)-SSC(I))*POIDS(I)                                             
      DO 260 J=1,5                                                              
      DO 265 K=J,5                                                              
      A(J,K)=A(J,K)+E(J)*E(K)                                                   
  265 CONTINUE                                                                  
  260 CONTINUE                                                                  
  230 CONTINUE                                                                  
      GOTO (270,275,280,285,290),IND                                            
  275 A(2,2)=1.                                                                 
      GOTO 285                                                                  
  270 A(2,2)=1.                                                                 
  280 A(3,3)=1.                                                                 
  285 A(4,4)=1.                                                                 
  290 DO 295 I=2,4                                                              
      L=I-1                                                                     
      DO 300 J=1,L                                                              
      A(I,J)=0.                                                                 
  300 CONTINUE                                                                  
  295 CONTINUE                                                                  
      A(1,1)=1./SQRT(A(1,1))                                                    
      DO 305 I=2,4                                                              
      A(1,I)=A(1,I)*A(1,1)                                                      
  305 CONTINUE                                                                  
      DO 310 I=2,4                                                              
      L=I-1                                                                     
      DO 315 J=1,L                                                              
      A(I,I)=A(I,I)-A(J,I)*A(J,I)                                               
      IF(A(I,I).LE.0.)GOTO580                                                   
  315 CONTINUE                                                                  
      RAIZDA=SQRT(A(I,I))                                                       
      A(I,I)=1./RAIZDA                                                          
      LL=I+1                                                                    
      DO 320 J=1,L                                                              
      DO 325 K=J,L                                                              
      A(I,J)=A(I,J)-A(K,I)*A(K,J)                                               
  325 CONTINUE                                                                  
      A(I,J)=A(I,J)*A(I,I)                                                      
  320 CONTINUE                                                                  
      IF (I.EQ.4) GOTO 310                                                      
      DO 330 J=LL,4                                                             
      DO 335 K=1,L                                                              
      A(I,J)=A(I,J)-A(K,I)*A(K,J)                                               
  335 CONTINUE                                                                  
      A(I,J)=A(I,J)*A(I,I)                                                      
  330 CONTINUE                                                                  
  310 CONTINUE                                                                  
      DO 337 J=1,4                                                              
      A(4,J)=A(4,J)*A(4,4)                                                      
  337 CONTINUE                                                                  
      DO 340 I=2,4                                                              
      IS=5-I                                                                    
      DIAG=A(IS,IS)                                                             
      DO 345 J=1,IS                                                             
      L=I-1                                                                     
      DO 350 K=1,L                                                              
      IM=5-K                                                                    
      A(IS,J)=A(IS,J)-A(IS,IM)*A(IM,J)                                          
  350 CONTINUE                                                                  
      A(IS,J)=A(IS,J)*DIAG                                                      
  345 CONTINUE                                                                  
  340 CONTINUE                                                                  
      DO 355 I=2,4                                                              
      L=I-1                                                                     
      DO 360 J=1,L                                                              
      A(J,I)=A(I,J)                                                             
  360 CONTINUE                                                                  
  355 CONTINUE                                                                  
      DO 365 I=1,4                                                              
      ECART(I)=0.                                                               
      DO 370 J=1,4                                                              
      ECART(I)=ECART(I)+A(I,J)*A(J,5)                                           
  370 CONTINUE                                                                  
  365 CONTINUE                                                                  
      GOTO (375,380,385,390,390),IND                                            
  375 ECART(2)=ECART(1)                                                         
  385 ECART(3)=ECART(1)                                                         
      GOTO 390                                                                  
  380 ECART(2)=ECART(1)                                                         
  390 DO 395 I=1,4                                                              
      PAR(I)=PAR(I)+ECART(I)                                                    
  395 CONTINUE                                                                  
      CALL PASAJE                                                               
      DO 400 J=1,4                                                              
      ECART(J)=REC(J)-DIR(J)                                                    
      DIR(J)=REC(J)                                                             
  400 CONTINUE                                                                  
      DO 410 J=1,4                                                              
      IF (ABS(ECART(J)).GT.ERROR(J)) GOTO 600                                   
  410 CONTINUE                                                                  
      DO 415 I=1,3                                                              
      PAR(I)=PAR(I)*BIDON1                                                      
      CALL PASAJE                                                               
      PAR(I)=PAR(I)*BIDON2                                                      
      DO 420 J=1,4                                                              
      B(I,J)=(REC(J)-DIR(J))*1000./PAR(I)                                       
  420 CONTINUE                                                                  
  415 CONTINUE                                                                  
      PAR(4)=PAR(4)+.01                                                         
      CALL PASAJE                                                               
      PAR(4)=PAR(4)-.01                                                         
      DO 425 J=1,4                                                              
      B(4,J)=(REC(J)-DIR(J))*100.                                               
  425 CONTINUE                                                                  
      JA=2                                                                      
      JB=2                                                                      
      GOTO (430,435,440,445,450),IND                                            
  440 JA=3                                                                      
  430 JB=3                                                                      
  435 DO 455 JV=JA,JB                                                           
      DO 460 J=1,4                                                              
      B(1,J)=B(1,J)+B(JV,J)                                                     
  460 CONTINUE                                                                  
  455 CONTINUE                                                                  
      DO 465 JV=JA,JB                                                           
      DO 470 J=1,4                                                              
      B(JV,J)=0.                                                                
  470 CONTINUE                                                                  
  465 CONTINUE                                                                  
  445 DO 475 J=1,4                                                              
      B(4,J)=0.                                                                 
  475 CONTINUE                                                                  
  450 RINDEP=A(5,5)/RINDEP                                                      
      DO 480 I=1,4                                                              
      ECART(I)=0.                                                               
      DO 485 J=1,4                                                              
      TEM(J)=0.                                                                 
      DO 490 K=1,4                                                              
      TEM(J)=TEM(J)+A(J,K)*B(K,I)                                               
  490 CONTINUE                                                                  
      ECART(I)=ECART(I)+TEM(J)*B(J,I)                                           
  485 CONTINUE                                                                  
      ECART(I)=SQRT(ECART(I)*RINDEP)                                            
  480 CONTINUE                                                                  
      GOTO (495,500,505,510,510),IND                                            
  495 ECART(2)=ECART(1)                                                         
  505 ECART(3)=ECART(2)                                                         
      ECART(2)=ECART(1)                                                         
      DIR(3)=DIR(2)                                                             
      DIR(2)=DIR(1)                                                             
      PAR(3)=PAR(2)                                                             
      PAR(2)=PAR(1)                                                             
      GOTO 510                                                                  
  500 ECART(2)=ECART(1)                                                         
C     I M P R E S S I O N   D E S   R E S U L T A T S                           
  510 X=PIRAD*DIR(4)                                                            
      SX=SIN(X)                                                                 
      KCLEF=1                                                                   
      Y=PIRAD*PAR(4)                                                            
      AR=PAR(1)*PAR(1)                                                          
      BR=PAR(2)*PAR(2)                                                          
      CR=PAR(3)*PAR(3)                                                          
      CAR=2.*PAR(3)*PAR(1)*COS(Y)                                               
  570 VOL=DIR(1)*DIR(2)*DIR(3)*SX                                               
      KLAVEQ=1                                                                  
      GOTO(571,572,573,574,575),IND                                             
  571 WRITE(6,1)DIR(1),VOL                                                      
      WRITE(85,1)DIR(1),VOL                                                    
      GOTO 576                                                                  
  572 WRITE(6,2)DIR(1),DIR(3),VOL                                               
      WRITE(85,2)DIR(1),DIR(3),VOL                                             
      GOTO 576                                                                  
  573 WRITE(6,3)DIR(1),DIR(3),VOL                                               
      WRITE(85,3)DIR(1),DIR(3),VOL                                             
      KLAVEQ=2                                                                  
      GOTO 576                                                                  
  574 WRITE(6,4)(DIR(I),I=1,3),VOL                                              
      WRITE(85,4)(DIR(I),I=1,3),VOL                                            
      GOTO 576                                                                  
  575 WRITE(6,5)(DIR(I),I=1,4),VOL                                              
      WRITE(85,5)(DIR(I),I=1,4),VOL                                            
      KLAVEQ=3                                                                  
  576 IF(KCLEF.NE.1) GOTO 579                                                   
      GOTO (578,577,577,578,578),IND                                            
  577 WRITE(6,6)ECART(1),ECART(3)                                               
      GOTO 579                                                                  
  578 WRITE(6,6)(ECART(I),I=1,NRIND)                                            
  579 IF(KDENS.EQ.1) GOTO 522                                                   
      WRITE(6,7)KZ                                                              
  522 WRITE(6,1004)                                                             
      DTHETA=0.                                                                 
      DELTAQ=0.                                                                 
      DO 525 I=1,N                                                              
      J=1                                                                       
      DIFQ=100.                                                                 
      KEY=1                                                                     
  555 XXH=XH(I,J)*XH(I,J)*AR                                                    
      XXK=XK(I,J)*XK(I,J)*BR                                                    
      XXL=XL(I,J)*XL(I,J)*CR                                                    
      CALQ=XXH+XXK+XXL                                                          
      GOTO (533,532,531),KLAVEQ                                                 
  531 CALQ=CALQ+(XH(I,J)*XL(I,J)*CAR)                                           
      GOTO 533                                                                  
  532 CALQ=CALQ+(XH(I,J)*XK(I,J)*AR)                                            
  533 DIFQ1=ABS(Q(I)-CALQ)                                                      
      DCAL=SQRT(CALQ)                                                           
      TCAL=2.*PIDEG*ASIN(WAVE2*DCAL)                                            
      DCAL=1./DCAL                                                              
      DIFD=D(I)-DCAL                                                            
      DIFT=TH(I)-TCAL                                                           
      GOTO (530,535),KEY                                                        
  535 WRITE(6,1005)XH(I,J),XK(I,J),XL(I,J),DCAL,DIFD,CALQ,TCAL,DIFT             
      GOTO 540                                                                  
  530 WRITE(6,1006)XH(I,J),XK(I,J),XL(I,J),D(I),DCAL,DIFD,EPSIL(I),             
     SQ(I),CALQ,TH(I),TCAL,DIFT                                                 
  540 IF (DIFQ1.GE.DIFQ) GOTO 545                                               
      DIFQ=DIFQ1                                                                
      DTET=ABS(DIFT)                                                            
  545 J=J+1                                                                     
      IF (QCAL(I,J).LT.0) GOTO 550                                              
      KEY=2                                                                     
      GOTO 555                                                                  
  550 DELTAQ=DIFQ+DELTAQ                                                        
      DTHETA=DTET+DTHETA                                                        
  525 CONTINUE                                                                  
C     F I G U R E S   D E   M E R I T E                                         
      DELTAQ=DELTAQ/N                                                           
      DTHETA=DTHETA/N                                                           
      FWOLFF=Q(N)/(2.*NPOSS*DELTAQ)                                             
      FINDEX=N/(NPOSS*DTHETA)
      WRITE(6,1007)N,NPOSS                                                      
      WRITE(6,1008)DELTAQ,DTHETA,DTH                                            
      WRITE(6,1009)N,FWOLFF,N,FINDEX,DTHETA,NPOSS                               
      WRITE(85,1007)N,NPOSS                                                 
      WRITE(85,1008)DELTAQ,DTHETA,DTH                                        
      WRITE(85,1009)N,FWOLFF,N,FINDEX,DTHETA,NPOSS                           
C ****************************************************************
C
C     LIMIT THE NUMBER OF ANSWERS FOR BATCH
C
C
      NUMANS=NUMANS+1
      IF(NUMANS.GT.10) THEN
                       WRITE(6,1267)
                       WRITE(85,1267)
1267  FORMAT(//' QUITE A FEW ANSWERS (10) HAVE BEEN FOUND !'/
     1   10X,'PLEASE EVALUATE THESE--THIS LIMIT COULD BE CHANGED.'/
     1   10X,'CALL JERRY JOHNSON FOR HELP.'//)
                        STOP 'ENDING JOB WITH 10 RESULTS'
      ENDIF
C
C
C
C *******************************************************************
      GOTO (560,580),LLAVE                                                      
CTEST N. 2 : FIGURES DE MERITE                                                  
  560 R=2.*DTH*NPOSS                                                            
      FWTEST=TAN(TH(N)/2.)/R                                                    
      FITEST=2.*N/R                                                             
      IF(FWTEST.GT.FWOLFF.OR.FITEST.GT.FINDEX) GOTO565                          
      GOTO 580                                                                  
C     SOLUTION DIFVOL                                                           
  565 WRITE(6,1010)                                                             
C             1010 : "SOLUTION NON AFFINEE"                                     
      LLAVE=2                                                                   
      KCLEF=2                                                                   
      DIR(1)=AA                                                                 
      DIR(2)=BB                                                                 
      DIR(3)=CC                                                                 
      DIR(4)=BETA                                                               
      X=PIRAD*BETA                                                              
      SX=SIN(X)                                                                 
      CX=COS(X)                                                                 
      SX2=SX*SX                                                                 
      AR=1./(AA*AA*SX2)                                                         
      BR=1./(BB*BB)                                                             
      CR=1./(CC*CC*SX2)                                                         
      CAR=-2.*CX/(AA*CC*SX2)                                                    
      GOTO 570                                                                  
    1 FORMAT(//20X,'S Y S T E M E   C U B I Q U E'                              
     S       //6X,'PARAMETRE DIRECT ',17X,'A=',F9.5,10X,'VOLUME=',F9.3)         
    2 FORMAT(//20X,'S Y S T E M E   T E T R A G O N A L'                        
     S       //6X,'PARAMETRES DIRECTS ',15X,'A=',F9.5,5X,'C=',F9.5,10X,         
     S       'VOLUME=',F9.3)                                                    
    3 FORMAT(//20X,'S Y S T E M E   H E X A G O N A L'                          
     S       //6X,'PARAMETRES DIRECTS ',15X,'A=',F9.5,5X,'C=',F9.5,10X,         
     S       'VOLUME=',F9.3/)                                                   
    4 FORMAT(//20X,'S Y S T E M E   O R T H O R H O M B I Q U E'                
     S       //6X,'PARAMETRES DIRECTS ',15X,'A=',F9.5,5X,'B=',F9.5,5X,          
     S       'C=',F9.5,10X,'VOLUME=',F9.3)                                      
    5 FORMAT(//20X,'S Y S T E M E   M O N O C L I N I Q U E'                    
     S       //6X,'PARAMETRES DIRECTS ',15X,'A=',F9.5,5X,'B=',F9.5,5X,          
     S       'C=',F9.5,5X,'BETA=',F7.3,10X,'VOLUME=',F9.3)                      
    6 FORMAT(6X,'ECARTS TYPES(S.D.)',18X,3(F9.5,7X),3X,F7.3/)                   
    7 FORMAT(/6X,'NOMBRE DE MOTIFS PAR MAILLE',7X,'Z=',I3/)                     
 1001 FORMAT(35X,3(F9.5,5X),F7.3,10X,F9.3,5X,'Z=',I2)                           
 1002 FORMAT(35X,3(F9.5,5X),F7.3,10X,F9.3)                                      
 1003 FORMAT(10X,'ECARTS-TYPES:',12X,3(F9.5,5X),F7.3)                           
 1004 FORMAT(/6X,'H',4X,'K',4X,'L',11X,'DOBS',6X,'DCAL',3X,'DOBS-DCAL'          
     S       ,2X,'EPSILON',8X,'QOBS',6X,'QCAL',6X,'2TH.OBS',2X,                 
     S       '2TH.CAL',3X,'DIF.2TH.'/)                                          
 1005 FORMAT(3X,3(F5.0),15X,2(F10.5),23X,F10.5,12X,2(F9.3))                     
 1006 FORMAT(3X,3(F5.0),5X,4(F10.5),3X,2(F10.5),3X,3(F9.3))                     
 1007 FORMAT(//10X,'* NOMBRE DE RAIES'/15X,'.- INTRODUITES =',I3/               
     S       15X,'.- CALCULEES =',I5)                                           
 1008 FORMAT(//10X,'* ECARTS MOYENS'/                                           
     S       25X,'<Q> =',E10.4/12X,'<DELTA(2-THETA)> =',E10.4/                  
     S       10X,'ERREUR MAX. TOLEREE=',E10.4)                                  
 1009 FORMAT(//10X,'* FIGURES DE MERITE'/17X,'1.- M(',I3,') =',F7.1,            
     S       25X,'(REF. 4)'/17X,'2.- F(',I3,') =',F7.1,'(',F6.4,',',I4,         
     S       ')',12X,'(REF. 5)'/)                                               
 1010 FORMAT(10X,'SOLUTION NON AFFINE'/)                                        
  580 V=VOL                                                                     
      RETURN                                                                    
      END                                                                       
      SUBROUTINE PASAJE                                                         
      COMMON D(50),Q(50),KQ(50),TH(50),EPSIL(50),EPSQ(50),KEPSQ(50),DTH,        
     1       V,VINF,VSUP,PAS,PAS2,PAS4,PAS8,PAS16,PAS32,PAS64,COEFF,RAP,        
     2       WAVE2,PIDEG,PIRAD,N,MH,MH2,MK,MK2,ML,ML2,KDENS,KZ,KQT(50)          
      COMMON QCAL(50,11),XH(50,11),XK(50,11),XL(50,11),AA,BB,CC,BETA            
      COMMON REC(4),PAR(4)                                                      
      SX=SIN(PAR(4)*PIRAD)                                                      
      VOL=1./(PAR(1)*PAR(2)*PAR(3)*SX)                                          
      REC(1)=PAR(2)*PAR(3)*VOL                                                  
      REC(2)=PAR(3)*PAR(1)*VOL*SX                                               
      REC(3)=PAR(1)*PAR(2)*VOL                                                  
      REC(4)=180.-PAR(4)                                                        
      RETURN                                                                    
      END                                                                       
