C***********************************************************************ðð00000
C*                                EVMI                                  0ð00000
C*                                                                      0ð00000
C*                HIROSHI TOKI             25.JUNE.1975                 0ð00000
C***********************************************************************ðð00000
C                                                                       0000000
C                                                                       0000000
C                                                                       0000000
C                    YADAV             JUNE 1976                        0000000
C          ************************************************             0000000
C*  1  16I5                                                             0ð00000
C*  IANUC  ATOMIC NUMBER                                                0ð00000
C*         IANUC > 0   VMI                                              0ð00000
C*         IANUC < 0   RIGID                                            0ð00000
C*         IANUC = 0   GO TO END                                        0ð00000
C*  IZNUC  CHARGE NUMBER                                                0ð00000
C*         IZNUC > 0   B(E2) AND Q-MOMENTS ARE CALCULATED               0ð00000
C*         IZNUC < 0   NOT CALCULATED                                   0ð00000
C*  IMAX   MAXIMUM SPIN   IMAX MUST BE SMALLER THAN 18.                 0ð00000
C*  NDATA  NUMBER OF DATA.   NDATA MUST BE SMALLER THAN 10.             0ð00000
C*  NALPH  NUMBER OF ALPHA.   NALPH MUST BE SMALLER THAN 10.            0ð00000
C*  IIMAX  MAXIMUM SPIN OF DATA USED FOR MINIMIZATION.                  0ð00000
C*  IOUT   CARD OUTPUT(2)  LIST OUTPUT(1)  NO(0)                        0000000
C*  2  8F10.5                                                           0ð00000
C*  GAMIP  GAMMA VALUE.   IF GAMIP < 0, GAM IS CALCULATED.              0ð00000
C#  BETA2                                                               0000000
C#  BETA4                                                               0000000
C#  BETA2 & BETA4 ARE Y2 AND Y4 DEFORMATIONS                            0000000
C*  3  8F10.5                                                           0ð00000
C*         (EEX(I),I=1,NDATA)                                           0ð00000
C*  EEX    EXPERIMENTAL EXCITATION ENERGIES. (MEV)                      0ð00000
C*  4  16(I3,I2)                                                        0ð00000
C*         ((IQUA(J,I),J=1,2),I=1,NDATA)                                0ð00000
C*  IQUA   IQUA(1,I)=SPIN                                               0ð00000
C*         IQUA(2,I)=ALPHA                                              0ð00000
C***********************************************************************ðð00000
C     EXTENDED VMI MODEL                                                0000000
C     DATE: 25,JUNE,1975                                                0000000
C     EVEN1: ASYMMETRIC CASE                                            0000000
C     DATE: 21,OCT,1974                                                 0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      EXTERNAL ECHS                                                     0000000
      COMMON/THET/THET1,THET2,THET3                                     0000000
      COMMON/CLECOM/BINOM(55,55)                                        0000000
      COMMON/ECHCOM/EEX(10),IQUA(2,10),NDATA,IIMAX                      0000000
      DIMENSION X(3),D(3),EN(10,18),EPS(10,18),EN1(10,18),EVEC(10,10,18)0000000
C     DIMENSION SUME(18)
      OPEN(5,FILE='triaxinp')                                            0000000
      OPEN(6,FILE='triaxout')
C                                                                       0000000
C     CALCULATION OF BINOMIAL FUNCTION                                  0000000
C                                                                       0000000
      DO 5 N1=1,55                                                      0000000
    5 BINOM(N1,N1)=1.D0                                                 0000000
      DO 6 N1=1,54                                                      0000000
      N2MAX=55-N1                                                       0000000
      DO 6 N2=1,N2MAX                                                   0000000
      N3=N1+N2                                                          0000000
      N4=N3-1                                                           0000000
    6 BINOM(N3,N1)=BINOM(N4,N1)*DFLOAT(N4)/DFLOAT(N2)                   0000000
C                                                                       0000000
C     INPUT                                                             0000000
C                                                                       0000000
   10 READ(5,20) IANUC,IZNUC,IMAX,NDATA,NALPH,IIMAX,IOUT                0000000
   20 FORMAT(16I5)                                                      0000000
      IF(IANUC.EQ.0) GO TO 1000                                         0000000
      READ(5,25) GAMIP,BET2,BET4                                        0000000
   25 FORMAT(8F10.5)                                                    0000000
      READ(5,30) (EEX(I),I=1,NDATA)                                     0000000
   30 FORMAT(8F10.5)                                                    0000000
      READ(5,35) ((IQUA(J,I),J=1,2),I=1,NDATA)                          0000000
   35 FORMAT(16(I3,I2))                                                 0000000
C                                                                       0000000
C     OUTPUT                                                            0000000
C                                                                       0000000
      WRITE(6,40) IANUC,IZNUC,IMAX,NDATA,NALPH,IIMAX,IOUT               0000000
   40 FORMAT('1IANUC=',I5,' IZNUC=',I5,'  IMAX=',I5,' NDATA=',I5,       0000000
     1  ' NALPH=',I5,' IIMAX=',I5,'  IOUT=',I5)                         0000000
      WRITE(6,45) GAMIP,BET2,BET4                                       0000000
   45 FORMAT(' GAMIP=',F10.5,' BETA2=',F10.5,' BETA4=',F10.5)           0000000
      WRITE(6,50) ((IQUA(J,I),J=1,2),I=1,NDATA)                         0000000
   50 FORMAT('0EXPERIMENTAL EXCITATION ENERGY (MEV)'/                   0000000
     1  1X,10(I5,I2,3X))                                                0000000
      WRITE(6,55) (EEX(I),I=1,NDATA)                                    0000000
   55 FORMAT(' ',10F10.5)                                               0000000
C     INITIAL VALUES                                                    0000000
      RAT2=EEX(2)/EEX(1)                                                0000000
      IF(RAT2.LT.2.D0) RAT2=2.D0                                        0000000
      SGAMA4=9.0/2.0*RAT2/(RAT2+1.0)**2                                 0000000
      SGAMA3=DSQRT(SGAMA4)                                              0000000
      GAMA3=DASIN(SGAMA3)                                               0000000
      GAMA=GAMA3/3.0                                                    0000000
      IF(GAMIP-0.00001) 70,70,75                                        0000000
   75 GAMA=GAMIP*3.141592654/180.0                                      0000000
      GAMA3=GAMA*3.0                                                    0000000
      SGAMA3=DSIN(GAMA3)                                                0000000
      SGAMA4=SGAMA3*SGAMA3                                              0000000
   70 CONTINUE                                                          0000000
      FK=10.0                                                           0000000
      EPS0=3.D0/EEX(1)                                                  0000000
      XPIE=3.141592654                                                  0000000
      BET22=BET2*BET2                                                   0000000
      BET42=BET4*BET4                                                   0000000
      BBET=5.0/4.0*BET42/BET22                                          0000000
      THET1=4.0/3.0*((DSIN(GAMA-2.0*XPIE/3.0))**2+BBET)                 0000000
      THET2=4.0/3.0*((DSIN(GAMA-4.0*XPIE/3.0))**2+BBET)                 0000000
      THET3=4.0/3.0*(DSIN(GAMA))**2                                     0000000
      WRITE(6,105)THET1,THET2,THET3,BBET,BET22,BET42                    0000000
  105 FORMAT('0 MOMENTS OF INERTIA *** THETA1=',F10.5,' THETA2=',F10.5, 0000000
     1 ' THETA3=',F10.5,' BETA4 TERM OF THETA1 & THETA2=',F10.5/        0000000
     2 ' SQUARE OF BETA2 AND BETA4 VALUES=',2F10.5)                     0000000
C***   IN ABOVE WE ASSUME  B4=0.5*B2, THEN I0=3B2*BET2**2               0000000
      IF(IANUC.LT.0) GO TO 100                                          0000000
      IF(EEX(2).LT.0.0) GO TO 510                                       0000000
C                                                                       0000000
C     SIMPLEX   (MINIMIZATION PROCEDURE)                                0000000
C                                                                       0000000
      X(1)=GAMA                                                         0000000
      X(2)=EPS0                                                         0000000
      X(3)=FK                                                           0000000
      D(1)=0.02                                                         0000000
      D(2)=2.0                                                          0000000
      D(3)=5.0                                                          0000000
      CALL SIMPLE(3,X,D,1.D0,0.5D0,2.D0,1.D-10,FECHS,IITE,ECHS)         0000000
      GAMA=X(1)                                                         0000000
      EPS0=X(2)                                                         0000000
      FK=X(3)                                                           0000000
      GAM=GAMA*180.0/3.141592654                                        0000000
      FC=1.0/FK                                                         0000000
      SIGMA=1.0/(2.0*FC*EPS0**3)                                        0000000
      WRITE(6,60) GAMA,EPS0,FK,GAM,FC,SIGMA,FECHS,IITE                  0000000
   60 FORMAT('0 GAMA=',F10.5,'  EPS0=',F10.5,'   FK=',F10.5,            0000000
     1  '   GAM=',F10.5,'    FC=',D10.4,' SIGMA=',D10.4,                0000000
     2  ' FECHS=',D10.4,'  IITE=',I5)                                   0000000
C                                                                       0000000
C     ENERGY WITH VMI                                                   0000000
C                                                                       0000000
      IANUC=-IANUC                                                      0000000
      CALL EVEN(GAMA,IMAX,EN,EVEC,IANUC)                                0000000
C     CALL VMI(IMAX,NALPH,GAMA,EPS0,FK,EN,EPS)                          0000000
      CALL VMIMOD(IMAX,NALPH,GAMA,EPS0,FK,EN,EPS)                       0000000
      GO TO 520                                                         0000000
  510 GAMA=GAMIP*3.141592654/180.0                                      0000000
      IANUC=-IANUC                                                      0000000
      CALL EVEN(GAMA,IMAX,EN,EVEC,IANUC)                                0000000
      CALL VMI3(EEX(1),EEX(3),EPS0,C,EN)                                0000000
      FK=1.0/C                                                          0000000
      SIGMA=1.0/(2.0*C*EPS0**3)                                         0000000
      WRITE(6,530) EPS0,C,SIGMA                                         0000000
  530 FORMAT('  EPS0=',F10.5,'     C=',F10.5,' SIGMA=',F10.5)           0000000
  520 WRITE(6,170) (J,J=1,NALPH)                                        0000000
  170 FORMAT('0   EXCITATION ENERGIES TAKING ACCOUNT VMI  (EPS)'        0000000
     1  //'  I ',I9,9I12)                                               0000000
      DO 180 I=2,IMAX                                                   0000000
      IEO=1-2*MOD(I,2)                                                  0000000
      IF(IEO) 192,1000,194                                              0000000
  194 ID=I/2+1                                                          0000000
      IALMX=MIN0(ID,NALPH)                                              0000000
      GO TO 196                                                         0000000
  192 ID=(I-1)/2                                                        0000000
      IALMX=MIN0(ID,NALPH-1)                                            0000000
  196 CONTINUE                                                          0000000
      IF(IALMX.EQ.0) GO TO 180                                          0000000
      DO 182 J=1,IALMX                                                  0000000
      IF(EEX(2).LT.0.0) EPS(J,I)=VMIEQ(EPS0,C,EN(J,I))                  0000000
      EN1(J,I)=EN(J,I)/(2.0*EPS(J,I))+(EPS(J,I)-EPS0)**2/(FK*2.D0)      0000000
  182 CONTINUE                                                          0000000
      WRITE(6,190) I,(EN1(J,I),J=1,IALMX)                               0000000
      WRITE(6,190) I,(EPS(J,I),J=1,IALMX)                               0000000
  190 FORMAT(' ',I2,3X,10F12.5)                                         0000000
  180 CONTINUE                                                          0000000
      GO TO 200                                                         0000000
C                                                                       0000000
C     ENERGY WITHOUT VMI                                                0000000
C                                                                       0000000
  100 CALL EVEN(GAMA,IMAX,EN,EVEC,IANUC)                                0000000
C     CHECK                                                             0000000
C     WRITE(6,101) EVEC(1,1,2),EVEC(1,2,2),EVEC(1,1,3),EVEC(1,1,4),     0000000
C    1  EVEC(1,2,4)                                                     0000000
C 101 FORMAT(' '5F10.5)                                                 0000000
C**   EN2=9.D0*(1.D0-DSQRT(1.D0-8.D0/9.D0*SGAMA4))/SGAMA4               0000000
C**   EPS(1,2)=3.D0/4.D0*EN2/EEX(1)                                     0000000
      ENX=1.0/THET1+1.0/THET2+1.0/THET3                                 0000000
      ENX2=ENX*ENX                                                      0000000
      ENY=(THET1+THET2+THET3)/(THET1*THET2*THET3)                       0000000
      EN2=ENX-DSQRT(ENX2-3.0*ENY)                                       0000000
      EPS(1,2)=EN2/EEX(1)                                               0000000
      GAM=GAMA*180.0/3.141592654                                        0000000
      ANUC=IANUC                                                        0000000
      BETA=408.*EPS(1,2)/((DABS(ANUC))**(7./3.))                        0000000
      WRITE(6,107)ENX,ENY,EN2,EEX(1)                                    0000000
  107 FORMAT('0 ENX,ENY,EN2,EEX(1) ARE=',4F10.5)                        0000000
      WRITE(6,110) GAMA,EPS(1,2),GAM,BETA                               0000000
  110 FORMAT('0 GAMA=',F10.5,'  EPS2=',F10.5,'   GAM=',F10.5,           0000000
     1  '  BETA=',F10.5)                                                0000000
      WRITE(6,270) (J,J=1,NALPH)                                        0000000
  270 FORMAT('0   EXCITATION ENERGIES NOT TAKING ACCOUNT VMI'//         0000000
     1'  I ',I9,9I12)                                                   0000000
      DO 280 I=2,IMAX                                                   0000000
      IEO=1-2*MOD(I,2)                                                  0000000
      IF(IEO) 284,1000,286                                              0000000
  286 ID=I/2+1                                                          0000000
      IALMX=MIN0(ID,NALPH)                                              0000000
      GO TO 288                                                         0000000
  284 ID=(I-1)/2                                                        0000000
      IALMX=MIN0(ID,NALPH-1)                                            0000000
  288 CONTINUE                                                          0000000
      DO 282 J=1,IALMX                                                  0000000
  282 EN1(J,I)=EN(J,I)/(2.0*EPS(1,2))                                   0000000
      WRITE(6,290) I,(EN1(J,I),J=1,IALMX)                               0000000
  290 FORMAT('0',I2,3X,10F12.5)                                         0000000
  280 CONTINUE                                                          0000000
      WRITE(6,291) (J,J=1,NALPH)                                        0000000
  291 FORMAT('0   EXCITATION ENERGIS IN UNIT OF 1/(2*EPS)'//            0000000
     1  '  I ',I9,9I12)                                                 0000000
      DO 292 I=2,IMAX                                                   0000000
      IEO=1-2*MOD(I,2)                                                  0000000
      IF(IEO) 294,1000,296                                              0000000
  296 ID=I/2+1                                                          0000000
      IALMX=MIN0(ID,NALPH)                                              0000000
      GO TO 298                                                         0000000
  294 ID=(I-1)/2                                                        0000000
      IALMX=MIN0(ID,NALPH-1)                                            0000000
  298 CONTINUE                                                          0000000
      WRITE(6,290) I,(EN(J,I),J=1,IALMX)                                0000000
C     SUME(I)=0.0                                                       0000000
C     DO 310 J=1,ID                                                     0000000
C 310 SUME(I)=SUME(I)+EN(J,I)                                           0000000
  292 CONTINUE                                                          0000000
C     WRITE(6,320) (I,I=2,IMAX)                                         0000000
C 320 FORMAT('0   SUMED VALUES OF THE EXCITATION ENERGIES'//            0000000
C    1  1X,10I12)                                                       0000000
C     WRITE(6,330) (SUME(I),I=2,IMAX)                                   0000000
C 330 FORMAT('0'5X,10F12.5)                                             0000000
  200 IANUC=IABS(IANUC)                                                 0000000
C                                                                       0000000
C     OUTPUT ENERGY AND WAVE FUNCTION ON CARDS                          0000000
C                                                                       0000000
      IF(IOUT.EQ.0) GO TO 410                                           0000000
      SQTWI=1.D0/DSQRT(2.D0)                                            0000000
      DO 420 J=1,2                                                      0000000
      IF(J.EQ.2.AND.EN1(J,2).GT.100.D0) GO TO 420                       0000000
      DO 4200 I=2,IMAX,2                                                0000000
      ID=I/2+1                                                          0000000
      DO 425 JJ=2,ID                                                    0000000
  425 EVEC(JJ,J,I)=EVEC(JJ,J,I)*SQTWI                                   0000000
      IF(IOUT.EQ.1) GO TO 426                                           0000000
      WRITE(7,430) EN1(J,I),(EVEC(JJ,J,I),JJ=1,ID)                      0000000
  426 WRITE(6,430) EN1(J,I),(EVEC(JJ,J,I),JJ=1,ID)                      0000000
 4200 CONTINUE 
  420 CONTINUE                                                          0000000
  430 FORMAT(3X,11F7.4)                                                 0000000
      IF(EN1(1,3).GT.100.D0) GO TO 410                                  0000000
      IMAX1=IMAX-1                                                      0000000
      DO 440 I=3,IMAX1,2                                                0000000
      ID=I/2                                                            0000000
      Y=0.0                                                             0000000
      DO 445 JJ=1,ID                                                    0000000
  445 EVEC(JJ,1,I)=EVEC(JJ,1,I)*SQTWI                                   0000000
      IF(IOUT.EQ.1) GO TO 440                                           0000000
      WRITE(7,430) EN1(1,I),Y,(EVEC(JJ,1,I),JJ=1,ID)                    0000000
  440 WRITE(6,430) EN1(1,I),Y,(EVEC(JJ,1,I),JJ=1,ID)                    0000000
  410 CONTINUE                                                          0000000
C                                                                       0000000
C     B(E2) VALUE AND Q-MOMENTS                                         0000000
C                                                                       0000000
      IF(IZNUC.LT.0) GO TO 10                                           0000000
      CALL BE2M(GAMA,BET2,IMAX,NALPH,EN1,EVEC,IANUC,IZNUC)              0000000
      CALL BE2T(GAMA,BET2,IMAX,NALPH,EN1,EVEC,IANUC,IZNUC)              0000000
      GO TO 10                                                          0000000
 1000 STOP                                                              0000000
      END                                                               0000000
      FUNCTION ECHS(X)                                                  0000000
C     CALCULATE THE KSI SQUARE VALUE                                    0000000
C     THIS SUBROUTINE IS THE EXTERNAL FUNCTION OF SIMPLE                0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      COMMON/ECHCOM/EEX(10),IQUA(2,10),NDATA,IIMAX                      0000000
      DIMENSION X(3),EN(10,18),EPS(10,18),ETH(10)                       0000000
      IANUC=1                                                           0000000
      GAMA=X(1)                                                         0000000
      EPS0=X(2)                                                         0000000
      FK=X(3)                                                           0000000
      GAM=GAMA*180.0/3.141592654                                        0000000
      IMAX=IIMAX                                                        0000000
      CALL EVEN(GAMA,IMAX,EN,EVEC,IANUC)                                0000000
      CALL VMIMOD(IMAX,2,GAMA,EPS0,FK,EN,EPS)                           0000000
C     CALL VMI(IMAX,2,GAMA,EPS0,FK,EN,EPS)                              0000000
      IF(EPS(1,2).LT.0.D0) GO TO 20                                     0000000
C     WRITE(6,70) EPS(2),EPS(4),EN(1,2),EN(2,2),EN(1,4)                 0000000
C  70 FORMAT('  EPS2=',F10.5,'  EPS4=',F10.5,'  EN21=',F10.5,           0000000
C    1  '  EN22=',F10.5,'  EN41=',F10.5)                                0000000
      DO 10 N=1,NDATA                                                   0000000
      J=IQUA(2,N)                                                       0000000
      I=IQUA(1,N)                                                       0000000
      ETH(N)=EN(J,I)/(EPS(J,I)*2.D0)+(EPS(J,I)-EPS0)**2/(FK*2.D0)       0000000
   10 CONTINUE                                                          0000000
      ECHS=0.0                                                          0000000
      DO 50 N=1,NDATA                                                   0000000
   50 ECHS=ECHS+(EEX(N)-ETH(N))**2                                      0000000
      EPS0MI=1.0/EEX(1)                                                 0000000
      EPS0MA=10.0*EPS0MI                                                0000000
      IF(EPS0.LT.EPS0MI) ECHS=ECHS+1.0-EPS0+EPS0MI                      0000000
      IF(EPS0.GT.EPS0MA) ECHS=ECHS+1.0+EPS0-EPS0MA                      0000000
      IF(GAMA.LT.0.D0) ECHS=ECHS+1.D0-GAMA                              0000000
      IF(GAMA.GT.30.D0) ECHS=ECHS+1.D0+GAMA-30.D0                       0000000
      IF(FK.LT.0.D0) ECHS=ECHS+FK*FK                                    0000000
C     WRITE(6,60) GAMA,EPS0,FK,GAM,ECHS,ETH(1),ETH(3)                   0000000
C  60 FORMAT('  GAMA=',F10.5,'  EPS0=',F10.5,'    FK=',F10.5,           0000000
C    1  '   GAM=',F10.5,'  ECHS=',F10.5,'  ETH1=',F10.5,                0000000
C    2  '  ETH3=',F10.5)                                                0000000
      RETURN                                                            0000000
   20 ECHS=10.0                                                         0000000
      RETURN                                                            0000000
      END                                                               0000000
      SUBROUTINE EVEN(GAMA,IMAX,EN,EVEC,IANUC)                          0000000
C     CALCULATE THE ENRGY AND THE WAVE FUNCTION OF EVEN SYSTEM          0000000
C     DAVYDOV AND FILLIPPOV                                             0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      COMMON/THET/THET1,THET2,THET3                                     0000000
      DIMENSION EN(10,18),EVEC(10,10,18)                                0000000
      DIMENSION FMAT(10,10),FMAT1(100),EIGV(100),EIG(10),EIG1(10)       0000000
      SSGAM3=DSIN(3.0*GAMA)                                             0000000
      SSGAM3=SSGAM3*SSGAM3                                              0000000
      SGAMA=DSIN(GAMA)                                                  0000000
      CGAMA=DCOS(GAMA)                                                  0000000
      SGAMA2=SGAMA*SGAMA                                                0000000
      GAMA1=3.0*(3.0-2.0*SGAMA2)/(3.0-4.0*SGAMA2)**2                    0000000
      GAMA2=-6.0*DSQRT(3.D0)*SGAMA*CGAMA/(3.0-4.0*SGAMA2)**2            0000000
      GAMA3=3.0/(4.0*SGAMA2)                                            0000000
      THETP=0.5*(1.0/THET1+1.0/THET2)                                   0000000
      THETM=0.5*(1.0/THET1-1.0/THET2)                                   0000000
      DO 50 J=1,10                                                      0000000
      DO 50 J1=1,10                                                     0000000
   50 FMAT(J,J1)=0.0                                                    0000000
C                                                                       0000000
C     BUILDING MATRIX ELEMENT                                           0000000
C                                                                       0000000
      DO 60 I=2,IMAX                                                    0000000
C     WRITE(6,65) I                                                     0000000
C  65 FORMAT('0I=',I5)                                                  0000000
      DO 52 J=1,10                                                      0000000
   52 EN(J,I)=0.0                                                       0000000
      IOE=1-2*MOD(I,2)                                                  0000000
      IF(IOE-1) 61,62,60                                                0000000
   61 ID=(I-1)/2                                                        0000000
      GO TO 63                                                          0000000
   62 ID=I/2+1                                                          0000000
   63 CONTINUE                                                          0000000
      FI=I                                                              0000000
      DO 70 J=1,ID                                                      0000000
      DO 70 J1=J,ID                                                     0000000
      IF(IOE+1) 81,82,81                                                0000000
   81 K=(J-1)*2                                                         0000000
      K1=(J1-1)*2                                                       0000000
      GO TO 83                                                          0000000
   82 K=J*2                                                             0000000
      K1=J1*2                                                           0000000
   83 FK=K                                                              0000000
      IF(K-K1) 80,90,80                                                 0000000
   80 IF(K+2-K1) 70,100,70                                              0000000
C* 90 FMAT(J,J1)=GAMA1*(FI*(FI+1.0)-FK*FK)+GAMA3*FK*FK                  0000000
   90 FMAT(J,J1)=THETP*(FI*(FI+1.0)-FK*FK)+GAMA3*FK*FK                  0000000
      GO TO 70                                                          0000000
C*100 FMAT(J,J1)=0.5*GAMA2*DSQRT((FI+FK+2.0)*(FI-FK-1.0)*(FI+FK+1.0)    0000000
C*   1*(FI-FK))                                                         0000000
  100 FMAT(J,J1)=0.5*THETM*DSQRT((FI+FK+2.0)*(FI-FK-1.0)*(FI+FK+1.0)    0000000
     1*(FI-FK))                                                         0000000
      IF(K.EQ.0) FMAT(J,J1)=DSQRT(2.D0)*FMAT(J,J1)                      0000000
   70 CONTINUE                                                          0000000
C                                                                       0000000
C     CHECK                                                             0000000
C                                                                       0000000
C     DO 72 J=1,ID                                                      0000000
C     WRITE(6,75) (FMAT(J,J1),J1=1,ID)                                  0000000
   75 FORMAT(' ',10E12.5)                                               0000000
   72 CONTINUE                                                          0000000
C                                                                       0000000
C     DIAGONALIZATION                                                   0000000
C                                                                       0000000
      K=0                                                               0000000
      DO 110 J1=1,ID                                                    0000000
      DO 110 J=1,ID                                                     0000000
      IF(J.GT.J1) GO TO 110                                             0000000
      K=K+1                                                             0000000
      FMAT1(K)=FMAT(J,J1)                                               0000000
  110 CONTINUE                                                          0000000
C                                                                       0000000
C      CHECK                                                            0000000
C                                                                       0000000
C     WRITE(6,115) (FMAT1(KK),KK=1,K)                                   0000000
  115 FORMAT(' FMAT1=',5D15.5)                                          0000000
      IKK=2                                                             0000000
      EIGV(1)=1.0                                                       0000000
      CALL MHS(FMAT1,ID,EIG,EIG1,IER)                                   0000000
      CALL MEVS(EIG,EIG1,ID,EIGV,ID,IER)                                0000000
      CALL MVS(FMAT1,ID,ID,EIGV,ID,IER)                                 0000000
      CALL ARANG(EIG,EIGV,ID,IKK)                                       0000000
C     WRITE(6,120) (EIG(J),J=1,ID)                                      0000000
C 120 FORMAT(' EIGENVALUE= ',5D15.5)                                    0000000
      DO 121 J=1,10                                                     0000000
  121 EN(J,I)=0.0                                                       0000000
      MID=MIN(10,ID)                                                    0000000
      DO 122 J=1,MID                                                    0000000
  122 EN(J,I)=EIG(J)                                                    0000000
      IF(IANUC.GT.0) GO TO 60                                           0000000
      DO 130 J=1,ID                                                     0000000
      DO 140 JJ=1,ID                                                    0000000
      JJ1=ID*(J-1)+JJ                                                   0000000
  140 EVEC(JJ,J,I)=EIGV(JJ1)                                            0000000
  130 CONTINUE                                                          0000000
   60 CONTINUE                                                          0000000
      RETURN                                                            0000000
      END                                                               0000000
      SUBROUTINE ARANG(A,R,N,IKK)                                       0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      DIMENSION A(1),R(1)                                               0000000
      IQ=-N                                                             0000000
      DO 1 I=1,N                                                        0000000
      IQ=IQ+N                                                           0000000
      JQ=N*(I-2)                                                        0000000
      DO 1 J=I,N                                                        0000000
      JQ=JQ+N                                                           0000000
      GO TO (5,6),IKK                                                   0000000
    5 IF(A(I)-A(J)) 2,1,1                                               0000000
    6 IF(A(I)-A(J)) 1,1,2                                               0000000
  2   X=A(I)                                                            0000000
      A(I)=A(J)                                                         0000000
      A(J)=X                                                            0000000
      DO 3 K=1,N                                                        0000000
      ILR=IQ+K                                                          0000000
      IMR=JQ+K                                                          0000000
      X=R(ILR)                                                          0000000
      R(ILR)=R(IMR)                                                     0000000
  3   R(IMR)=X                                                          0000000
  1   CONTINUE                                                          0000000
      RETURN                                                            0000000
      END                                                               0000000
C     *** MHS ********* VERSION 1, MODIFICATION LEVEL 0 *** DKO10331 ***0000000
C     *                                                                *0000000
C     *   REDUCE A SYMMETRIC MATRIX TO SYMMETRIC TRIDIAGONAL FORM      *0000000
C     *                                                                *0000000
C     *   5736-XM7 COPYRIGHT IBM CORP. 1971                            *0000000
C     *   REFER TO INSTRUCTIONS ON COPYRIGHT NOTICE FORM NO. 120-2083  *0000000
C     *   FE SERVICE NO. 200281                                        *0000000
C     *                                                                *0000000
C     ******************************************************************0000000
C                                                                       0000000
      SUBROUTINE MHS(A,N,D,CD,IER)                                      0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      DIMENSION A(100),D(10),CD(10)                                     0000000
      CALL COIBM                                                        0000000
      ISW  =IER                                                         0000000
      IER  =0                                                           0000000
      IF (N) 1,1,3                                                      0000000
    1 IER  =1000                                                        0000000
      IF (ISW+12345) 2,23,2                                             0000000
    2 CALL WIER(IER,10331)                                              0000000
      GOTO 23                                                           0000000
    3 D(1) =A(1)                                                        0000000
      IF (N-3) 4,6,6                                                    0000000
    4 IF (N-2) 23,5,23                                                  0000000
    5 D(2) =A(3)                                                        0000000
      CD(2)=A(2)                                                        0000000
      GOTO 23                                                           0000000
    6 E1   =1.D-12                                                      0000000
      ID   =1                                                           0000000
      ICD  =2                                                           0000000
      LU   =N-1                                                         0000000
      DO 22 M=2,LU                                                      0000000
         S    =0.0D0                                                    0000000
C                                                                       0000000
C     NORM OF THE PART OF THE M-TH ROW WHICH IS TO BE REDUCED           0000000
C                                                                       0000000
         J    =ICD+M                                                    0000000
         LL   =M+1                                                      0000000
         DO 7 I=LL,N                                                    0000000
            D(I) =A(J)                                                  0000000
            DT   =A(J)                                                  0000000
            S    =S+DT*DT                                               0000000
    7       J    =J+I                                                   0000000
         DT   =A(ICD)                                                   0000000
         DT   =DT*DT                                                    0000000
         IF (S-DT*E1) 8,8,9                                             0000000
C                                                                       0000000
C     BYPASS M-TH SIMILARITY                                            0000000
C                                                                       0000000
    8    CD(M)=A(ICD)                                                   0000000
         A(ID)=0.0D0                                                    0000000
         GOTO 21                                                        0000000
C                                                                       0000000
C     NEW CODIAGONAL ELEMENT                                            0000000
C                                                                       0000000
    9    DT   =DSQRT(S+DT)                                              0000000
         IF (A(ICD)) 11,11,10                                           0000000
   10    DT   =-DT                                                      0000000
   11    A(ID)=DT                                                       0000000
         D(M) =A(ICD)-A(ID)                                             0000000
         J    =ICD                                                      0000000
         DT   =0.0D0                                                    0000000
C                                                                       0000000
C     COMPUTE VECTORS DEFINING THE TRANSFORMATION                       0000000
C                                                                       0000000
         DO 16 L=M,N                                                    0000000
            S    =0.0D0                                                 0000000
            LK   =J                                                     0000000
            DO 12 K=M,L                                                 0000000
               LK   =LK+1                                               0000000
               DD   =D(K)                                               0000000
   12          S    =S+A(LK)*DD                                         0000000
            LK   =LK+L                                                  0000000
            LL   =L+1                                                   0000000
            IF (N-LL) 15,13,13                                          0000000
   13       DO 14 K=LL,N                                                0000000
               DD   =D(K)                                               0000000
               S    =S+A(LK)*DD                                         0000000
   14          LK   =LK+K                                               0000000
   15       DT   =DT+S*D(L)                                             0000000
            CD(L)=S                                                     0000000
   16       J    =J+L                                                   0000000
         DT   =5.D-1*DT                                                 0000000
         T    =D(M)*A(ID)                                               0000000
         DO 17 L=M,N                                                    0000000
            D(L) =D(L)/T                                                0000000
   17       CD(L)=CD(L)+DT*D(L)                                         0000000
         J    =ICD                                                      0000000
C                                                                       0000000
C     PERFORM M-TH SIMILARITY                                           0000000
C                                                                       0000000
         DO 20 K=M,N                                                    0000000
            LK   =J                                                     0000000
            IF (K-M) 20,18,18                                           0000000
   18       DO 19 L=M,K                                                 0000000
               LK   =LK+1                                               0000000
               S    =A(LK)                                              0000000
               DD   =CD(K)                                              0000000
               DB   =CD(L)                                              0000000
               S    =S+D(L)*DD+D(K)*DB                                  0000000
   19          A(LK)=S                                                  0000000
   20       J    =J+K                                                   0000000
         CD(M)=A(ID)                                                    0000000
   21    ID   =ICD+1                                                    0000000
         D(M) =A(ID)                                                    0000000
   22    ICD  =ID+M                                                     0000000
      CD(N)=A(ICD)                                                      0000000
      A(ID)=A(ICD)                                                      0000000
      D(N) =A(ICD+1)                                                    0000000
      CD(1)=0.0D0                                                       0000000
   23 RETURN                                                            0000000
      END                                                               0000000
C     *** MEVS ******** VERSION 1, MODIFICATION LEVEL 0 *** DKO10332 ***0000000
C     *                                                                *0000000
C     *   EIGENVALUES AND EIGENVECTORS OF A                            *0000000
C     *   REAL SYMMETRIC TRIDIAGONAL MATRIX                            *0000000
C     *   QR-METHOD - IMPLICIT SHIFT                                   *0000000
C     *                                                                *0000000
C     *   5736-XM7 COPYRIGHT IBM CORP. 1971                            *0000000
C     *   REFER TO INSTRUCTIONS ON COPYRIGHT NOTICE FORM NO. 120-2083  *0000000
C     *   FE SERVICE NO. 200281                                        *0000000
C     *                                                                *0000000
C     ******************************************************************0000000
C                                                                       0000000
      SUBROUTINE MEVS(D,CD,N,EV,IEV,IER)                                0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      DIMENSION D(1),CD(1),EV(1)                                        0000000
      EQUIVALENCE (SEV,ISEV)                                            0000000
      CALL COIBM                                                        0000000
      ISW  =IER                                                         0000000
      IER  =0                                                           0000000
      IF (N) 1,1,2                                                      0000000
    1 IER  =1000                                                        0000000
      GOTO 33                                                           0000000
    2 SEV  =EV(1)                                                       0000000
      IF (ISEV) 3,12,3                                                  0000000
    3 IF (IABS(IEV)-N) 1,4,4                                            0000000
    4 INEV =1                                                           0000000
      IF (IEV) 5,6,6                                                    0000000
    5 IREV =-IEV                                                        0000000
      ICEV =INEV                                                        0000000
      IER  =1                                                           0000000
      GOTO 7                                                            0000000
    6 IREV =INEV                                                        0000000
      ICEV =IEV                                                         0000000
    7 IV1  =1                                                           0000000
      DO 11 I=1,N                                                       0000000
         IV   =IV1                                                      0000000
         DO 10 J=1,N                                                    0000000
            IF (I-J) 8,9,8                                              0000000
    8       EV(IV)=0.D0                                                 0000000
            GOTO 10                                                     0000000
    9       EV(IV)=1.D0                                                 0000000
   10       IV   =IV+IREV                                               0000000
   11    IV1  =IV1+ICEV                                                 0000000
   12 E1   =1.D-6                                                       0000000
      MIT  =20                                                          0000000
      M    =N                                                           0000000
   13 M1   =M-1                                                         0000000
      IF (M1) 32,32,14                                                  0000000
C                                                                       0000000
C     START LOOP FOR ITERATION                                          0000000
C                                                                       0000000
   14 DO 30 IT=1,MIT                                                    0000000
         U    =DABS(D(M))                                               0000000
         J    =M                                                        0000000
         DO 15 JJ=2,M                                                   0000000
            T    =U                                                     0000000
            J1   =J-1                                                   0000000
            U    =DABS(D(J1))                                           0000000
            IF (DABS(CD(J))-E1*(T+U)) 16,16,15                          0000000
   15       J    =J1                                                    0000000
   16    IF (J-M) 17,31,17                                              0000000
   17    CDJ  =CD(J)                                                    0000000
         S    =CD(M)                                                    0000000
C                                                                       0000000
C     DETERMINE THE SHIFT                                               0000000
C                                                                       0000000
         T    =(D(M)-D(M1))*.5                                          0000000
         U    =DSQRT(T*T+S*S)                                           0000000
         IF (T) 18,19,19                                                0000000
   18    U    =-U                                                       0000000
   19    Q    =D(J)-(U+T+D(M1))                                         0000000
C                                                                       0000000
C     QR-TRANSFORMATION                                                 0000000
C                                                                       0000000
         C    =1.D0                                                     0000000
         S    =1.D0                                                     0000000
         U    =0.D0                                                     0000000
         I1   =J                                                        0000000
         I2   =I1+1                                                     0000000
         IF (M-I2) 29,20,20                                             0000000
   20    IF (ISEV) 21,22,21                                             0000000
   21    INIT =I1*ICEV+1                                                0000000
   22    DO 28 I=I2,M                                                   0000000
            H    =C*CD(I)                                               0000000
            S    =S*CD(I)                                               0000000
            IF (DABS(S)-DABS(Q)) 23,24,24                               0000000
   23       S    =S/Q                                                   0000000
            R    =DSQRT(S*S+1.D0)                                       0000000
            CD(I1)=Q*R                                                  0000000
            C    =1.D0/R                                                0000000
            S    =S*C                                                   0000000
            GOTO 25                                                     0000000
   24       C    =Q/S                                                   0000000
            R    =DSQRT(C*C+1.D0)                                       0000000
            CD(I1)=S*R                                                  0000000
            S    =1.D0/R                                                0000000
            C    =C*S                                                   0000000
   25       G    =D(I1)-U                                               0000000
            T    =(D(I)-G)*S+H*(C+C)                                    0000000
            U    =S*T                                                   0000000
            D(I1)=G+U                                                   0000000
            Q    =C*T-H                                                 0000000
            IF (ISEV) 26,28,26                                          0000000
C                                                                       0000000
C     FORM EIGENVECTORS                                                 0000000
C                                                                       0000000
   26       IV   =INIT                                                  0000000
            DO 27 K=1,N                                                 0000000
               IV1  =IV-ICEV                                            0000000
               T    =EV(IV1)                                            0000000
               R    =EV(IV)                                             0000000
               EV(IV1)=C*T+S*R                                          0000000
               EV(IV)=C*R-S*T                                           0000000
   27          IV   =IV+IREV                                            0000000
            INIT =INIT+ICEV                                             0000000
   28       I1   =I                                                     0000000
C                                                                       0000000
C     END QR-TRANSFORMATION                                             0000000
C                                                                       0000000
   29    CD(J)=CDJ                                                      0000000
         CD(M)=Q                                                        0000000
   30    D(M) =D(M)-U                                                   0000000
C                                                                       0000000
C     END LOOP FOR ITERATION                                            0000000
C                                                                       0000000
      IER  =IER+2000                                                    0000000
      CD(1)=N-M                                                         0000000
      GOTO 33                                                           0000000
   31 M    =M-1                                                         0000000
      GOTO 13                                                           0000000
   32 IF (IER) 33,35,33                                                 0000000
   33 IF (ISW+12345) 34,35,34                                           0000000
   34 CALL WIER(IER,10332)                                              0000000
   35 RETURN                                                            0000000
      END                                                               0000000
C     *** COIBM ******* VERSION 1, MODIFICATION LEVEL 0 *** DKO00000 ***0000000
C     *                                                                *0000000
C     *   COPYRIGHT STATEMENT FOR PROGRAM PRODUCT SL-MATH              *0000000
C     *                                                                *0000000
C     *   5736-XM7 COPYRIGHT IBM CORP. 1971                            *0000000
C     *   REFER TO INSTRUCTIONS ON COPYRIGHT NOTICE FORM NO. 120-2083  *0000000
C     *   FE SERVICE NO. 200281                                        *0000000
C     *                                                                *0000000
C     ******************************************************************0000000
C                                                                       0000000
      SUBROUTINE COIBM                                                  0000000
    1 FORMAT ('5736-XM7 COPYRIGHT IBM CORP. 1971')                      0000000
      RETURN                                                            0000000
      END                                                               0000000
C     *** WIER ******** VERSION 1, MODIFICATION LEVEL 0 *** DKO00001 ***0000000
C     *                                                                *0000000
C     *   WRITE ERROR MESSAGES FOR SL-MATH SUBROUTINES                 *0000000
C     *                                                                *0000000
C     *   5736-XM7 COPYRIGHT IBM CORP. 1971                            *0000000
C     *   REFER TO INSTRUCTIONS ON COPYRIGHT NOTICE FORM NO. 120-2083  *0000000
C     *   FE SERVICE NO. 200281                                        *0000000
C     *                                                                *0000000
C     ******************************************************************0000000
C                                                                       0000000
      SUBROUTINE WIER(IER,NO)                                           0000000
      CALL COIBM                                                        0000000
      WRITE (6,1) NO,IER                                                0000000
    1 FORMAT (//' ***** DKO',I5,' RAISED ERROR INDICATOR TO ',I4,       0000000
     1     ' *****'///)                                                 0000000
      RETURN                                                            0000000
      END                                                               0000000
C     *** MVS ********* VERSION 1, MODIFICATION LEVEL 0 *** DKO10336 ***0000000
C     *                                                                *0000000
C     *   BACKTRANSFORMATION OF THE EIGENVECTORS OF A                  *0000000
C     *   REAL SYMMETRIC MATRIX                                        *0000000
C     *                                                                *0000000
C     *   5736-XM7 COPYRIGHT IBM CORP. 1971                            *0000000
C     *   REFER TO INSTRUCTIONS ON COPYRIGHT NOTICE FORM NO. 120-2083  *0000000
C     *   FE SERVICE NO. 200281                                        *0000000
C     *                                                                *0000000
C     ******************************************************************0000000
C                                                                       0000000
      SUBROUTINE MVS(A,N,M,EV,IEV,IER)                                  0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      DIMENSION A(1),EV(1)                                              0000000
      CALL COIBM                                                        0000000
      ISW  =IER                                                         0000000
      IER  =0                                                           0000000
      INEV =1                                                           0000000
      IF (N) 1,1,2                                                      0000000
    1 IER  =1000                                                        0000000
      GOTO 18                                                           0000000
    2 IF (M) 1,1,3                                                      0000000
    3 IF (N-M) 1,4,4                                                    0000000
    4 IF (IABS(IEV)-N) 1,5,5                                            0000000
    5 IF (IEV) 6,7,7                                                    0000000
    6 IREV =-IEV                                                        0000000
      ICEV =INEV                                                        0000000
      IER  =1                                                           0000000
      GOTO 8                                                            0000000
    7 IREV =INEV                                                        0000000
      ICEV =IEV                                                         0000000
    8 N1   =N-1                                                         0000000
      N2   =N-2                                                         0000000
      IV1  =1                                                           0000000
      IDN1 =N*N1/2                                                      0000000
      DO 17 L=1,M                                                       0000000
         IF (N2) 14,14,9                                                0000000
    9    ID   =IDN1                                                     0000000
         KK   =N2*IREV+IV1                                              0000000
         DO 13 K=1,N2                                                   0000000
            NMK  =N-K                                                   0000000
            ICD  =ID-1                                                  0000000
            ID   =ID-NMK                                                0000000
            IF (A(ID)) 10,13,10                                         0000000
C                                                                       0000000
C     ORTHOGONAL TRANSFORMATIONS                                        0000000
C                                                                       0000000
   10       C    =A(ICD)-A(ID)                                          0000000
            S    =0.D0                                                  0000000
            J    =ICD                                                   0000000
            IV   =KK                                                    0000000
            DO 11 I=NMK,N                                               0000000
               DP   =A(J)                                               0000000
               S    =S+DP*EV(IV)                                        0000000
               J    =J+I                                                0000000
   11          IV   =IV+IREV                                            0000000
            S    =S/A(ID)                                               0000000
            T    =(S-EV(KK))/C                                          0000000
            EV(KK)=S                                                    0000000
            J    =ICD+NMK                                               0000000
            K1   =NMK+1                                                 0000000
            IV   =KK+IREV                                               0000000
            DO 12 I=K1,N                                                0000000
               EV(IV)=EV(IV)+T*A(J)                                     0000000
               J    =J+I                                                0000000
   12          IV   =IV+IREV                                            0000000
   13       KK   =KK-IREV                                               0000000
   14    S    =0.D0                                                     0000000
         IV   =IV1                                                      0000000
C                                                                       0000000
C     NORMALIZATION                                                     0000000
C                                                                       0000000
         DO 15 I=1,N                                                    0000000
            DP   =EV(IV)                                                0000000
            S    =DP*DP+S                                               0000000
   15       IV   =IV+IREV                                               0000000
         C    =DSQRT(S)                                                 0000000
         IV   =IV1                                                      0000000
         DO 16 I=1,N                                                    0000000
            EV(IV)=EV(IV)/C                                             0000000
   16       IV   =IV+IREV                                               0000000
   17    IV1  =IV1+ICEV                                                 0000000
      IF (IER) 18,20,18                                                 0000000
   18 IF (ISW+12345) 19,20,19                                           0000000
   19 CALL WIER(IER,10336)                                              0000000
   20 RETURN                                                            0000000
      END                                                               0000000
      SUBROUTINE VMI(IMAX,NALPH,GAMA,EPS0,FK,EN,EPS)                    0000000
C     EVMI                                                              0000000
C     PROGRAMED BY TOKI                                                 0000000
C     DATE: 11,11,1974                                                  0000000
C     VMI IN ASYMMETRIC ROTOR                                           0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      DIMENSION EPS(10,18),EN(10,18)                                    0000000
C     ***  STEP5  VMI  ***                                              0000000
      EPMIN=EPS0*2.D0/3.D0                                              0000000
      EPMAX=EPS0*50                                                     0000000
      DO 120 I=2,IMAX                                                   0000000
      ID=NALPH                                                          0000000
      IEO=1-2*MOD(I,2)                                                  0000000
      IF(IEO.EQ.-1) ID=NALPH-1                                          0000000
      DO 120 J=1,ID                                                     0000000
      IF(EN(J,I).LT.0.00001D0) GO TO 120                                0000000
      DELEP=(EPMAX-EPMIN)/10.0                                          0000000
      EPS1=EPMIN                                                        0000000
      K=0                                                               0000000
  110 K=K+1                                                             0000000
      FEPS=EPS1**3-EPS0*EPS1*EPS1-EN(J,I)*FK/2.D0                       0000000
      IF(DABS(FEPS).LT.1.D-10) GO TO 80                                 0000000
      IF(DELEP.LT.1.D-10) GO TO 80                                      0000000
      IF(FEPS) 90,80,100                                                0000000
   90 EPS1=EPS1+DELEP                                                   0000000
      IF(EPS1.GT.EPMAX) GO TO 1010                                      0000000
      GO TO 110                                                         0000000
  100 DELEP=DELEP/2.0                                                   0000000
      EPS1=EPS1-DELEP                                                   0000000
      GO TO 110                                                         0000000
   80 EPS(J,I)=EPS1                                                     0000000
  120 CONTINUE                                                          0000000
C     WRITE(6,130) (EPS(I),I=2,IMAX)                                    0000000
C 130 FORMAT('0EPS=',5D15.7)                                            0000000
      RETURN                                                            0000000
 1010 WRITE(6,1020)  EPMIN,EPMAX,EPS1,EPS0,FK,EN(J,I)                   0000000
 1020 FORMAT('0THE SOLUTION EXCEED OVER THE MAXIMAM VALUE.'             0000000
     1  //' EPMIN=',F10.5,' EPMAX=',F10.5,'  EPS1=',F10.5,              0000000
     2  '  EPS0=',F10.5,'    FK=',F10.5,'    EN=',F10.5)                0000000
      EPS(1,2)=-1.0                                                     0000000
      RETURN                                                            0000000
      END                                                               0000000
      SUBROUTINE VMIMOD(IMAX,NALPH,GAMA,EPS0,FK,EN,EPS)                 0000000
C     DATE : 23. 10. 1975                                               0000000
C     VMI IN ASYMMETRIC ROTOR                                           0000000
C     KARDANISCHE FORMEL (MATHEMATISCHE FORMELSAMMLUNG BY KARL          0000000
C     HOTTMANN    P.15)                                                 0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      DIMENSION EPS(10,18),EN(10,18)                                    0000000
C     ***  STEP5  VMI  ***                                              0000000
      B=-EPS0                                                           0000000
      P=-B*B/9.D0                                                       0000000
      P3=P*P*P                                                          0000000
      DO 120 I=2,IMAX                                                   0000000
      ID=NALPH                                                          0000000
      IEO=1-2*MOD(I,2)                                                  0000000
      IF(IEO.EQ.-1) ID=NALPH-1                                          0000000
      DO 120 J=1,ID                                                     0000000
      IF(EN(J,I).LT.0.00001D0) GO TO 120                                0000000
      D=-EN(J,I)*FK/2.D0                                                0000000
      Q=B*B*B/27.D0+D/2.D0                                              0000000
      PPQQ=Q*Q+P3                                                       0000000
      IF(PPQQ.LT.0.D0) GO TO 1010                                       0000000
      SQPQ=DSQRT(PPQQ)                                                  0000000
      U=DCBRT(-Q+SQPQ)                                                  0000000
      V=DCBRT(-Q-SQPQ)                                                  0000000
      Y=U+V                                                             0000000
      EPS(J,I)=Y-B/3.D0                                                 0000000
  120 CONTINUE                                                          0000000
      RETURN                                                            0000000
 1010 WRITE(6,1020) PPQQ                                                0000000
 1020 FORMAT('  THERE ARE THREE REAL SOLUTIONS.  PPQQ=',F10.5)          0000000
      EPS(1,2)=-1.0                                                     0000000
      RETURN                                                            0000000
      END                                                               0000000
      SUBROUTINE BE2M(GAMA,BET2,IMAX,NALPH,EN1,EVEC,IANUC,IZNUC)        0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      COMMON/CLECOM/BINOM(55,55)                                        0000000
      DIMENSION EVEC(10,10,18),EN1(10,18)                               0000000
C     DIMENSION SUMQ(18)                                                0000000
C     CHECK                                                             0000000
C     WRITE(6,101) EVEC(1,1,2),EVEC(1,2,2),EVEC(1,1,3),EVEC(1,1,4),     0000000
C    1  EVEC(1,2,4)                                                     0000000
C 101 FORMAT(' '5F10.5)                                                 0000000
      SGAMA=DSIN(GAMA)                                                  0000000
      CGAMA=DCOS(GAMA)                                                  0000000
      ZNUC=IZNUC                                                        0000000
      ANUC=IANUC                                                        0000000
      PAI=3.141592654                                                   0000000
      GAM=GAMA*180.0/PAI                                                0000000
      Q0=3.D0*ZNUC*1.2*1.2*ANUC**(2.0/3.0)/DSQRT(PAI*5.D0)*BET2         0000000
      A0=DSQRT(5.D0/(16.0*PAI))                                         0000000
      WRITE(6,10) GAM,Q0                                                0000000
   10 FORMAT('1  E2 PROPERTIES OF ASYMMETRIC ROTOR',5X,                 0000000
     1  'SPECTROSCOPIC QUADRUPOLE MOMENTS'//'   GAM=',F10.5,            0000000
     2  '    Q0=',F10.5//'  I IAL   EN       RE2',10X,'SQM',10X,        0000000
     3  'SQM0'/19X,'E*BARN',6X,'E*FM**2'/)                              0000000
      DO 20 I=2,IMAX                                                    0000000
C     SUMQ(I)=0.0                                                       0000000
      IEO=1-2*MOD(I,2)                                                  0000000
      IF(IEO) 30,1000,40                                                0000000
   40 ID=I/2+1                                                          0000000
      IALMX=MIN0(ID,NALPH)                                              0000000
      GO TO 50                                                          0000000
   30 ID=(I-1)/2                                                        0000000
      NALPH1=NALPH-1                                                    0000000
      IALMX=MIN0(ID,NALPH1)                                             0000000
   50 DO 20 IAL=1,IALMX                                                 0000000
      II=2*I                                                            0000000
      RE2=0.0                                                           0000000
      DO 60 N=1,ID                                                      0000000
      IF(IEO) 70,1000,80                                                0000000
   80 IO=2*(N-1)                                                        0000000
      GO TO 90                                                          0000000
   70 IO=2*N                                                            0000000
   90 IIO=2*IO                                                          0000000
      IIIO=-IIO                                                         0000000
      DO 60 N1=1,ID                                                     0000000
      IF(IEO) 100,1000,110                                              0000000
  110 IO1=2*(N1-1)                                                      0000000
      GO TO 120                                                         0000000
  100 IO1=2*N1                                                          0000000
  120 IIO1=2*IO1                                                        0000000
      IIIO1=-IIO1                                                       0000000
      CALL CLEBSH(CG1,II,4,II,IIO1,0,IIO)                               0000000
      CALL CLEBSH(CG2,II,4,II,IIO1,4,IIO)                               0000000
      CALL CLEBSH(CG3,II,4,II,IIIO1,4,IIO)                              0000000
      CALL CLEBSH(CG4,II,4,II,IIIO1,4,IIIO)                             0000000
      SCG=CGAMA*CG1+SGAMA/DSQRT(2.D0)*(CG2+CG3*DFLOAT(IEO)+CG4)         0000000
      IF(IO) 1000,130,140                                               0000000
  130 IF(IO1.NE.0) SCG=SCG/DSQRT(2.D0)                                  0000000
      GO TO 60                                                          0000000
  140 IF(IO1.EQ.0) SCG=SCG/DSQRT(2.D0)                                  0000000
   60 RE2=RE2+EVEC(N,IAL,I)*EVEC(N1,IAL,I)*SCG                          0000000
      RE2=RE2*Q0*A0                                                     0000000
C     RE2=RE2*SQRT(DFLOAT(2*I+1))/100.
      RE2=RE2*DSQRT(DFLOAT(2*I+1))/100.
      CALL CLEBSH(CG,II,4,II,II,0,II)                                   0000000
      SQM=CG*RE2/A0                                                     0000000
      SQM0=SQM/Q0                                                       0000000
C     SUMQ(I)=SUMQ(I)+SQM0                                              0000000
      WRITE(6,160) I,IAL,EN1(IAL,I),RE2,SQM,SQM0                        0000000
      WRITE(9,190) I,IAL,I,IAL,RE2                                      0000000
  160 FORMAT(' ',2I3,F10.5,3D13.5)                                      0000000
  190 FORMAT('23.   ',4I5,F10.5)                                        0000000
   20 CONTINUE                                                          0000000
C     WRITE(6,170) (I,I=2,IMAX)                                         0000000
C 170 FORMAT('0   SUMED VALUES OF THE QUADRUPOLE MOMENTS'//             0000000
C    1  1X,10I12)                                                       0000000
C     WRITE(6,180) (SUMQ(I),I=2,IMAX)                                   0000000
C 180 FORMAT('0',5X,10F12.5)                                            0000000
      RETURN                                                            0000000
 1000 STOP                                                              0000000
      END                                                               0000000
      SUBROUTINE BE2T(GAMA,BET2,IMAX,NALPH,EN1,EVEC,IANUC,IZNUC)        0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      COMMON/CLECOM/BINOM(55,55)                                        0000000
      DIMENSION EVEC(10,10,18),EN1(10,18)                               0000000
C     CHECK                                                             0000000
C     WRITE(6,101) EVEC(1,1,2),EVEC(1,2,2),EVEC(1,1,3),EVEC(1,1,4),     0000000
C    1  EVEC(1,2,4)                                                     0000000
C 101 FORMAT(' ',5F10.5)                                                0000000
      SGAMA=DSIN(GAMA)                                                  0000000
      CGAMA=DCOS(GAMA)                                                  0000000
      ZNUC=IZNUC                                                        0000000
      ANUC=IANUC                                                        0000000
      PAI=3.141592654                                                   0000000
      GAM=GAMA*180.0/PAI                                                0000000
      Q0=3.D0*ZNUC*1.2*1.2*ANUC**(2.0/3.0)/DSQRT(PAI*5.D0)*BET2         0000000
      A0=DSQRT(5.D0/(16.0*PAI))                                         0000000
      TUNT=1.23D9                                                       0000000
      WRITE(6,10) GAM,Q0                                                0000000
   10 FORMAT('1  E2 PROPERTIES OF ASYMMETRIC ROTOR',5X,                 0000000
     1  'ELECTRIC QUADRUPOLE TRANSITIONS'//'   GAM=',F10.5,             0000000
     2  '    Q0=',F10.5//'  I IAL   EN',5X,'  I IAL   EN',6X,           0000000
     3  'RE2',10X,'BE2',10X,'TE2',10X,'BE20'/5X,'INITIAL',10X,          0000000
     4  'FINAL',8X,'E*BARN',6X,'(E*FM**2)**2',1X,'1/SEC'/)              0000000
      DO 20 I=2,IMAX                                                    0000000
      IEO=1-2*MOD(I,2)                                                  0000000
      IF(IEO) 30,1000,40                                                0000000
   40 ID=I/2+1                                                          0000000
      IALMX=MIN0(ID,NALPH)                                              0000000
      GO TO 50                                                          0000000
   30 ID=(I-1)/2                                                        0000000
      NALPH1=NALPH-1                                                    0000000
      IALMX=MIN0(ID,NALPH1)                                             0000000
   50 DO 20 IAL=1,IALMX                                                 0000000
      II=2*I                                                            0000000
      DO 20 IN=1,3                                                      0000000
      I1=I-3+IN                                                         0000000
      IF(I1.EQ.1) GO TO 20                                              0000000
      I1EO=1-2*MOD(I1,2)                                                0000000
      IF(I1EO) 35,1000,45                                               0000000
   45 ID1=I1/2+1                                                        0000000
      IALMX1=MIN0(ID1,NALPH)                                            0000000
      IF(I1.EQ.I) IALMX1=IAL-1                                          0000000
      GO TO 55                                                          0000000
   35 ID1=(I1-1)/2                                                      0000000
      NALPH1=NALPH-1                                                    0000000
      IALMX1=MIN0(ID1,NALPH1)                                           0000000
      IF(I1.EQ.I) IALMX1=IAL-1                                          0000000
   55 IF(IALMX1.EQ.0) GO TO 20                                          0000000
      DO 25 IAL1=1,IALMX1                                               0000000
      II1=2*I1                                                          0000000
      RE2=0.0                                                           0000000
      DO 60 N=1,ID                                                      0000000
      IF(IEO) 70,1000,80                                                0000000
   80 IO=2*(N-1)                                                        0000000
      GO TO 90                                                          0000000
   70 IO=2*N                                                            0000000
   90 IIO=2*IO                                                          0000000
      IIIO=-IIO                                                         0000000
      DO 60 N1=1,ID1                                                    0000000
      IF(I1EO) 100,1000,110                                             0000000
  110 IO1=2*(N1-1)                                                      0000000
      GO TO 120                                                         0000000
  100 IO1=2*N1                                                          0000000
  120 IIO1=2*IO1                                                        0000000
      IIIO1=-IIO1                                                       0000000
      CALL CLEBSH(CG1,II1,4,II,IIO1,0,IIO)                              0000000
      CALL CLEBSH(CG2,II1,4,II,IIO1,4,IIO)                              0000000
      CALL CLEBSH(CG3,II1,4,II,IIIO1,4,IIO)                             0000000
      CALL CLEBSH(CG4,II1,4,II,IIIO1,4,IIIO)                            0000000
      SCG=CGAMA*CG1+SGAMA/DSQRT(2.D0)*(CG2+CG3*DFLOAT(I1EO)+CG4*        0000000
     1  DFLOAT(IEO*I1EO))                                               0000000
      IF(IO) 1000,130,140                                               0000000
  130 IF(IO1.NE.0) SCG=SCG/DSQRT(2.D0)                                  0000000
      GO TO 150                                                         0000000
  140 IF(IO1.EQ.0) SCG=SCG/DSQRT(2.D0)                                  0000000
  150 IF(I1) 151,152,151                                                0000000
  151 RE2=RE2+EVEC(N,IAL,I)*EVEC(N1,IAL1,I1)*SCG                        0000000
      GO TO 60                                                          0000000
  152 RE2=RE2+EVEC(N,IAL,I)*SCG                                         0000000
   60 CONTINUE                                                          0000000
      RE2=RE2*Q0*A0                                                     0000000
      BE2=RE2*RE2*DFLOAT(2*I1+1)/DFLOAT(2*I+1)                          0000000
C     RE2=RE2*SQRT(DFLOAT(2*I1+1))/100.
      RE2=RE2*DSQRT(DFLOAT(2*I1+1))/100.
      RE2=RE2*(-1.)**(I-I1)
      IF(I1) 155,156,155                                                0000000
  155 TE2=TUNT*DABS(EN1(IAL,I)-EN1(IAL1,I1))**5*BE2                     0000000
      GO TO 157                                                         0000000
  156 TE2=TUNT*EN1(IAL,I)**5*BE2                                        0000000
  157 BE20=BE2/(Q0*Q0)*16.0*PAI                                         0000000
      IF(I1.EQ.0) GO TO 165                                             0000000
      WRITE(6,160) I1,IAL1,EN1(IAL1,I1),I,IAL,EN1(IAL,I),               0000000
     1  RE2,BE2,TE2,BE20                                                0000000
      WRITE(9,180) I1,IAL1,I,IAL,RE2                                    0000000
      GO TO 25                                                          0000000
  165 ZERO=0.0                                                          0000000
      WRITE(6,160) I1,IAL1,ZERO,I,IAL,EN1(IAL,I),                       0000000
     1  RE2,BE2,TE2,BE20                                                0000000
      WRITE(9,180) I1,IAL1,I,IAL,RE2                                    0000000
  160 FORMAT(' ',2I3,F10.5,2I3,F10.5,4D13.5)                            0000000
  180 FORMAT('23.    ',4I5,F10.5)                                       0000000
   25 CONTINUE                                                          0000000
   20 CONTINUE                                                          0000000
      RETURN                                                            0000000
 1000 STOP                                                              0000000
      END                                                               0000000
      SUBROUTINE SIMPLE(N,X,D,ALPHA,BETA,GAMMA,EPS,FNEU,JZAHL,FKT)      0000000
      IMPLICIT REAL *8 (A-H,O-Z)                                        0000001
      DIMENSION X(3),XX(3,4),YY(3),XNULL(3),XR(3),XE(3),XC(3)           0000002
      DIMENSION F(3),D(3)                                               0000003
      NP1=N+1                                                           0000004
      ALPH1=1+ALPHA                                                     0000005
      BETA1=1-BETA                                                      0000006
      GAMM1=1-GAMMA                                                     0000007
C                                                                       0000008
C     ERSTELLEN DER KOEFFIZIENTENMATRIX                                 0000009
C                                                                       0000010
      JZAHL=0                                                           0000011
      DO 1 K=1,NP1                                                      0000012
      DO 1 I=1,N                                                        0000013
      XX(I,K)=X(I)                                                      0000014
      IF(K.EQ.I) XX(I,K)=XX(I,K)+D(I)                                   0000015
    1 CONTINUE                                                          0000016
C                                                                       0000017
C     BERECHNUNG DER DAZUGEHOERIGEN FUNKTIONSWERTEN                     0000018
C                                                                       0000019
      DO 3 I=1,NP1                                                      0000020
      DO 4 K=1,N                                                        0000021
    4 YY(K)=XX(K,I)                                                     0000022
    3 F(I)=FKT(YY)                                                      0000023
C                                                                       0000024
C     MINIMUM                                                           0000025
   55 LL=1                                                              0000026
      FL=F(1)                                                           0000027
      DO 5 I=2,NP1                                                      0000028
      IF((F(I)-FL).GE.0.) GO TO 5                                       0000029
      FL=F(I)                                                           0000030
      LL=I                                                              0000031
    5 CONTINUE                                                          0000032
C                                                                       0000033
C     MAXIMUM                                                           0000034
   66 LH=1                                                              0000035
      FH=F(1)                                                           0000036
      DO 6 I=2,NP1                                                      0000037
      IF((F(I)-FH).LE.0.) GO TO 6                                       0000038
      FH=F(I)                                                           0000039
      LH=I                                                              0000040
    6 CONTINUE                                                          0000041
C                                                                       0000042
C     ZWEITES MAXIMUM                                                   0000043
   77 II=2                                                              0000044
      LH2=1                                                             0000045
      FH2=F(1)                                                          0000046
      IF(LH-1) 74,75,74                                                 0000047
   75 LH2=2                                                             0000048
      II=3                                                              0000049
      FH2=F(2)                                                          0000050
   74 DO 7 I=II,NP1                                                     0000051
      IF(I.EQ.LH) GO TO 7                                               0000052
      IF((F(I)-FH2).LE.0.) GO TO 7                                      0000053
      FH2=F(I)                                                          0000054
      LH2=I                                                             0000055
    7 CONTINUE                                                          0000056
      JZAHL=JZAHL+1                                                     0000057
      FQUER=0.                                                          0000058
      DO 8 K=1,NP1                                                      0000059
    8 FQUER=FQUER+F(K)                                                  0000060
      FQUER=FQUER/NP1                                                   0000061
      S=0.                                                              0000062
      DO 9 I=1,NP1                                                      0000063
    9 S=S+(F(I)-FQUER)**2                                               0000064
      S=DSQRT(S/N)                                                      0000065
      IF(S.LE.EPS)   GO TO 216                                          0000066
      GO TO 213                                                         0000067
216   DO 11 I=1,N                                                       0000068
11    X(I)=XX(I,LL)                                                     0000069
      FNEU=F(LL)                                                        0000070
       RETURN                                                           0000071
C                                                                       0000072
C                                                                       0000073
  213 DO 46 I=1,N                                                       0000074
      XNULL(I)=0.                                                       0000075
      DO 47 K=1,NP1                                                     0000076
      IF(K.EQ.LH) GO TO 47                                              0000077
      XNULL(I)=XNULL(I)+XX(I,K)                                         0000078
   47 CONTINUE                                                          0000079
   46 XNULL(I)=XNULL(I)/N                                               0000080
      DO 13 K=1,N                                                       0000081
   13 XR(K)=ALPH1*XNULL(K)-ALPHA*XX(K,LH)                               0000082
      FR=FKT(XR)                                                        0000083
      IF(FR.GE.F(LL)) GO TO 26                                          0000084
      DO 16 K=1,N                                                       0000085
   16 XE(K)=GAMMA*XR(K)+GAMM1*XNULL(K)                                  0000086
      FE=FKT(XE)                                                        0000087
      IF(FE.GE.F(LL)) GO TO 23                                          0000088
      DO 19 I=1,N                                                       0000089
   19 XX(I,LH)=XE(I)                                                    0000090
      F(LH)=FE                                                          0000091
   20 LL=LH                                                             0000092
      LH=LH2                                                            0000093
      GO TO 77                                                          0000094
   23 DO 24 I=1,N                                                       0000095
   24 XX(I,LH)=XR(I)                                                    0000096
      F(LH)=FR                                                          0000097
      GO TO 20                                                          0000098
   26 IF(FR.GT.F(LH2)) GO TO 30                                         0000099
      DO 27 I=1,N                                                       0000100
   27 XX(I,LH)=XR(I)                                                    0000101
      F(LH)=FR                                                          0000102
      LH=LH2                                                            0000103
      GO TO 77                                                          0000104
   30 IF(FR.GE.F(LH)) GO TO 32                                          0000105
      DO 31 I=1,N                                                       0000106
   31 XX(I,LH)=XR(I)                                                    0000107
      F(LH)=FR                                                          0000108
   32 DO 33 I=1,N                                                       0000109
   33 XC(I)=BETA*XX(I,LH)+BETA1*XNULL(I)                                0000110
      FC=FKT(XC)                                                        0000111
      IF(FC.GE.F(LH)) GO TO 36                                          0000112
      DO 35 I=1,N                                                       0000113
   35 XX(I,LH)=XC(I)                                                    0000114
      F(LH)=FC                                                          0000115
      GO TO 55                                                          0000116
   36 DO 42 I=1,N                                                       0000117
      DO 42 K=1,NP1                                                     0000118
      IF(K.EQ.LL) GO TO  42                                             0000119
       XX(I,K)=(XX(I,K)+XX(I,LL))*0.5                                   0000120
   42 CONTINUE                                                          0000121
      DO 37 I=1,NP1                                                     0000122
      DO 38 K=1,N                                                       0000123
   38 YY(K)=XX(K,I)                                                     0000124
   37 F(I)=FKT(YY)                                                      0000125
      GO TO  66                                                         0000126
      END                                                               0000127
      SUBROUTINE CLEBSH(CG,JA,JB,JC,MA,MB,MC)                           0000000
C     ALL THE ARGUMENT IS TWICE OF THE ANGULAR MOMRNTUM.                0000000
C     CHECKED ON 26.10.1974.                                            0000000
C     THE CG IS EQUAL TO BRINK SATCHLER.                                0000000
C     P34 (2.34)                                                        0000000
C                                                                       0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      COMMON/CLECOM/BINOM(55,55)                                        0000000
      CG=0.0                                                            0000000
      MAB=MA+MB                                                         0000000
      IF(MAB.NE.MC) GO TO 20                                            0000000
      IF(JA+JB.LT.JC) GO TO 20                                          0000000
      IF(JB+JC.LT.JA) GO TO 20                                          0000000
      IF(JC+JA.LT.JB) GO TO 20                                          0000000
      J=(JA+JB+JC)/2                                                    0004410
      J1=J-JA                                                           0004420
      J2=J-JB                                                           0004430
      J3=J-JC                                                           0004440
      M1=(JA-MA)/2                                                      0004450
      M2=(JB-MB)/2                                                      0004460
      M3=(JC-MC)/2                                                      0004470
      N1=M1+MA                                                          0004480
      N2=M2+MB                                                          0004490
      N3=M3+MC                                                          0004500
      CG=0.0                                                            0004510
      KMIN=MAX0(0,M1-J2,N2-J1)                                          0004520
      KMAX=MIN0(J3,M1,N2)                                               0004530
      IF(KMIN-KMAX) 12,12,20                                            0004540
   12 K1=J-N3                                                           0004550
      K2=J-M3                                                           0004560
      TT=BINOM(J+1,N3+1)*BINOM(K1+1,M1+1)                               0004570
      TT=BINOM(J+1,M3+1)*BINOM(K2+1,N1+1)*TT                            0004580
      TT=BINOM(J+1,JC+1)*BINOM(JC+1,J1+1)/TT                            0004590
      CEF=DFLOAT(JC+1)*TT/DFLOAT(J+1)                                   0004600
      CEF=DSQRT(CEF)                                                    0004610
      S=(-1.0)**KMIN                                                    0004620
      K=KMIN                                                            0004630
   15 L1=M1-K                                                           0004640
      L2=N2-K                                                           0004650
      CG=CG+S*BINOM(J3+1,K+1)*BINOM(J2+1,L1+1)*BINOM(J1+1,L2+1)         0004660
      S=-S                                                              0004670
      K=K+1                                                             0004680
      IF(K.LE.KMAX) GO TO 15                                            0004690
      CG=CEF*CG                                                         0004700
   20 RETURN                                                            0004710
      END                                                               0004720
      FUNCTION VMIEQ(EPS0,C,EN)                                         0000000
C     DATE : 23. 10. 1975                                               0000000
C     VMIEQ IS MADE AT 13. NOV. 1975                                    0000000
C     VMI IN ASYMMETRIC ROTOR                                           0000000
C     KARDANISCHE FORMEL (MATHEMATISCHE FORMELSAMMLUNG BY KARL          0000000
C     ROTTMANN    P.15) AND CASUS IRREDUCIBILIS                         0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      IF(C) 50,50,40                                                    0000000
   50 IF(EN.LT.-EPS0**3*8.0*C/27.0) GO TO 40                            0000000
      VMIEQ=EPS0*2.0/3.0                                                0000000
      RETURN                                                            0000000
   40 FK=1.0/C                                                          0000000
      B=-EPS0                                                           0000000
      P=-B*B/9.D0                                                       0000000
      P3=P*P*P                                                          0000000
      D=-EN*FK/2.D0                                                     0000000
      Q=B*B*B/27.D0+D/2.D0                                              0000000
      PPQQ=Q*Q+P3                                                       0000000
      IF(PPQQ) 10,20,20                                                 0000000
C     P**3+Q**2 IS POSITIVE                                             0000000
   20 SQPQ=DSQRT(PPQQ)                                                  0000000
      U=DCBRT(-Q+SQPQ)                                                  0000000
      V=DCBRT(-Q-SQPQ)                                                  0000000
      Y=U+V                                                             0000000
      VMIEQ=Y-B/3.D0                                                    0000000
C     WRITE(6,30) VMIEQ,EPS0,C,EN                                       0000000
C  30 FORMAT(' VMIEQ',10F10.5)                                          0000000
      RETURN                                                            0000000
C     P**3+Q**2 IS NEGATIVE                                             0000000
   10 COSPHI=-Q/DSQRT(-P3)                                              0000000
      PHI=DACOS(COSPHI)                                                 0000000
      PHI3=PHI/3.D0                                                     0000000
      PAI=3.141592654D0                                                 0000000
      DP=DSQRT(-P)                                                      0000000
      Y1=2.D0*DP*DCOS(PHI3)                                             0000000
      Y2=-2.D0*DP*DCOS(PHI3+PAI/3.D0)                                   0000000
      Y3=-2.D0*DP*DCOS(PHI3-PAI/3.D0)                                   0000000
      Y=DMAX1(Y1,Y2,Y3)                                                 0000000
      VMIEQ=Y-B/3.D0                                                    0000000
C     WRITE(6,30) VMIEQ,EPS0,C,EN                                       0000000
      RETURN                                                            0000000
      END                                                               0000000
      FUNCTION DCBRT(X)                                                 0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      Y=DABS(X)                                                         0000000
      IF(Y.LT.0.0000000001) GO TO 1                                     0000000
      PH=X/Y                                                            0000000
      X=Y                                                               0000000
      DCBRT=X**(1.0/3.0)                                                0000000
      DCBRT=DCBRT*PH                                                    0000000
      RETURN                                                            0000000
    1 DCBRT=0.0                                                         0000000
      RETURN                                                            0000000
      END                                                               0000000
      SUBROUTINE VMI3(EN21,EN41,EPS0,C,EN)                              0000000
C     EVMI                                                              0000000
C     PROGRAMED BY TOKI                                                 0000000
C     DATE: 11,11,1974                                                  0000000
C     VMI IN ASYMMETRIC ROTOR                                           0000000
      IMPLICIT REAL*8 (A-H,O-Z)                                         0000000
      DIMENSION EN(10,18)                                               0000000
      DELTA=1.D-10                                                      0000000
      EPS2=EN(1,2)/(2.0*EN21)                                           0000000
      EPSLIM=EPS2*2.0                                                   0000000
      EPS4=EN(1,4)/(2.0*EN41)                                           0000000
      WRITE(6,13) EN(1,2),EN(1,4)                                       0000000
   13 FORMAT('0  EN21=',F10.5,'   EN41=',F10.5)                         0000000
C     ***  STEP  EPS0,C,SIGM  ***                                       0000000
      EPS0=(10.0*EPS2**3-3.0*EPS4**3)/(10.0*EPS2*EPS2-3.0*EPS4*EPS4)    0000000
      C=3.0/(EPS2**3-EPS0*EPS2*EPS2)                                    0000000
      KK=0                                                              0000000
  110 EPS2=VMIEQ(EPS0,C,EN(1,2))                                        0000000
      EPS4=VMIEQ(EPS0,C,EN(1,4))                                        0000000
      WRITE(6,11) EPS0,C,EPS2,EPS4                                      0000000
   11 FORMAT('  EPS0=',F10.5,' C=',F10.5,' EPS2=',F10.5,' EPS4=',F10.5) 0000000
      EPSDIF=DABS(3.D0*EPS2-2.D0*EPS0)                                  0000000
      IF(EPSDIF-DELTA) 30,30,40                                         0000000
   30 KK=KK+1                                                           0000000
      EPS0=EPS0+0.5                                                     0000000
      C=C+0.005                                                         0000000
      IF(KK.GT.30) GO TO 140                                            0000000
      GO TO 110                                                         0000000
   40 DEE2=EPS2/(3.0*EPS2-2.0*EPS0)                                     0000000
      DEE4=EPS4/(3.0*EPS4-2.0*EPS0)                                     0000000
      DCE2=EN(1,2)/(2.0*C*C)/(3.0*EPS2*EPS2-2.0*EPS2*EPS0)              0000000
      DCE4=EN(1,4)/(2.0*C*C)/(3.0*EPS4*EPS4-2.0*EPS4*EPS0)              0000000
      F=EN(1,2)/(2.0*EPS2)+(EPS2-EPS0)**2*C/2.0-EN21                    0000000
      G=EN(1,4)/(2.0*EPS4)+(EPS4-EPS0)**2*C/2.0-EN41                    0000000
      FE=-EN(1,2)*DEE2/(2.0*EPS2*EPS2)+C*(EPS2-EPS0)*(DEE2-1.0)         0000000
      GE=-EN(1,4)*DEE4/(2.0*EPS4*EPS4)+C*(EPS4-EPS0)*(DEE4-1.0)         0000000
      FC=-EN(1,2)*DCE2/(2.0*EPS2*EPS2)+(EPS2-EPS0)**2/2.0               0000000
     1  +C*(EPS2-EPS0)*DCE2                                             0000000
      GC=-EN(1,4)*DCE4/(2.0*EPS4*EPS4)+(EPS4-EPS0)**2/2.0               0000000
     1  +C*(EPS4-EPS0)*DCE4                                             0000000
      DENOM=FE*GC-GE*FC                                                 0000000
      DELE=(G*FC-F*GC)/DENOM                                            0000000
      DELC=(F*GE-G*FE)/DENOM                                            0000000
      WRITE(6,12) DELE,DELC,F,G                                         0000000
   12 FORMAT(' DELE=',F10.5,' DELC=',F10.5,' F=',F10.5,' G=',F10.5)     0000000
      IF(DABS(F).LT.DELTA.AND.DABS(G).LT.DELTA) GO TO 120               0000000
      EPS0=EPS0+DELE                                                    0000000
      C=C+DELC                                                          0000000
      WRITE(6,10) EPS0,C,EPS2,EPS4                                      0000000
   10 FORMAT(' EPS0=',F10.5,' C=',F10.5,' EPS2=',F10.5,' EPS4=',F10.5)  0000000
      IF(DABS(EPS0)-EPSLIM) 110,130,130                                 0000000
  130 EPS0=EPS0-DELE+0.5                                                0000000
      C=C-DELC+0.005                                                    0000000
      KK=KK+1                                                           0000000
      IF(KK.GT.30) GO TO 140                                            0000000
      GO TO 110                                                         0000000
  120 RETURN                                                            0000000
  140 WRITE(6,150) EPS0,C                                               0000000
  150 FORMAT('0 VMI3 CAN NOT FIND THE SOLUTION'/5X,' EPS0=',F10.5,      0000000
     1  ' C=',F10.5)                                                    0000000
      EPS0=EPS0-2.0                                                     0000000
      IF(EPS0.GT.-EPSLIM) GO TO 110                                     0000000
      STOP                                                              0000000
      END                                                               0000000
