engylvl.bas

Newman程序是一套很好的光谱软件计算程序

REM    PROGRAM ENGYLVL
REM    TO CONSTRUCT THE ENERGY MATRIX FROM CF PARAMETERS
REM      AND FIND ITS EIGENVALUES AND EIGENVECTORS
REM
REM    "CRYSTAL FIELD HANDBOOK" (SECTION 3.1.3 AND APPENDIX 2)
REM    EDITED BY D.J. NEWMAN AND BETTY NG
REM    CAMBRIDGE UNIVERSITY PRESS

       DEFINT I-J, M-N
       DEFSNG A-E, G-H, S-T, V, Z
       DEFDBL R
       DIM B2Q(-2 TO 2), B4Q(-4 TO 4), B6Q(-6 TO 6)
       DIM RMC(1 TO 3)
       DIM EMAT(1 TO 17, 1 TO 17), A(1 TO 17, 1 TO 17)
       DIM D(1 TO 17), V(1 TO 17, 1 TO 17), B(1 TO 100), Z(1 TO 100)
DECLARE SUB threej (R0#, AJ1!, AJ2!, AJ3!, AM1!, AM2!, AM3!)

       PRINT "THIS PROGRAM CALCULATES THE ENERGY LEVELS OF A J-MULTIPLET"
       PRINT "  FROM A SET OF CRYSTAL FIELD PARAMETERS"
       PRINT "   "
       INPUT "J VALUE (ANY REAL NO. <=8) OF THE MULTIPLET = "; AJ1
       PRINT "PROVIDE REDUCED MATRIX ELEMENTS FOR THIS MULTIPLET"
       INPUT "REDUCED MATRIX ELEMENT FOR RANK 2 = "; RMC(1)
       INPUT "REDUCED MATRIX ELEMENT FOR RANK 4 = "; RMC(2)
       INPUT "REDUCED MATRIX ELEMENT FOR RANK 6 = "; RMC(3)
       PRINT "    "

      INPUT "FILENAME OF CRYSTAL FIELD PARAMETERS = "; FILE1$
      PRINT "   "
      OPEN FILE1$ FOR INPUT AS #5

      INPUT #5, DUMMY$
      INPUT #5, B2Q(-2), B2Q(-1), B2Q(0), B2Q(1), B2Q(2)
      INPUT #5, DUMMY$
      INPUT #5, B4Q(-4), B4Q(-3), B4Q(-2), B4Q(-1), B4Q(0), B4Q(1), B4Q(2), B4Q(3), B4Q(4)
      INPUT #5, DUMMY$
      INPUT #5, B6Q(-6), B6Q(-5), B6Q(-4), B6Q(-3), B6Q(-2), B6Q(-1), B6Q(0), B6Q(1), B6Q(2), B6Q(3), B6Q(4), B6Q(5), B6Q(6)
      CLOSE #5

REM   INITIALISE THE ENERGY MATRIX
      FOR I = 1 TO 17
          FOR J = 1 TO 17
              EMAT(I, J) = 0!
          NEXT J
      NEXT I
REM   CONSTRUCT THE ENERGY MATRIX
      PRINT "CONSTRUCTION OF ENERGY MATRIX BEGINS"
      PRINT "PLEASE BE PATIENT. IT MAY TAKE A WHILE ....."
      FOR AM1 = -AJ1 TO AJ1 STEP 1
          IA1 = AM1 + AJ1 + 1
          FOR AM2 = -AJ1 TO AJ1 STEP 1
              IA2 = AM2 + AJ1 + 1
REM   SUMMING OVER RANK 2
              FOR M3 = 0 TO 2
                 AM3 = M3
                  IF M3 = 0 THEN
                   CALL threej(R0, AJ1, 2!, AJ1, -AM1, AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1)) * R0 * RMC(1) * B2Q(M3)
                   ELSE
                   CALL threej(R0, AJ1, 2!, AJ1, -AM1, AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1)) * R0 * RMC(1) * B2Q(M3)
                   CALL threej(R0, AJ1, 2!, AJ1, -AM1, -AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1 + M3)) * R0 * RMC(1) * B2Q(M3)
                  END IF
              NEXT M3
REM  SUMMING OVER RANK 4
              FOR M3 = 0 TO 4
                AM3 = M3
                  IF M3 = 0 THEN
                   CALL threej(R0, AJ1, 4!, AJ1, -AM1, AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1)) * R0 * RMC(2) * B4Q(M3)
                   ELSE
                   CALL threej(R0, AJ1, 4!, AJ1, -AM1, AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1)) * R0 * RMC(2) * B4Q(M3)
                   CALL threej(R0, AJ1, 4, AJ1, -AM1, -AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1 + M3)) * R0 * RMC(2) * B4Q(M3)
                  END IF
              NEXT M3
REM  SUMMING OVER RANK 6
              FOR M3 = 0 TO 6
              AM3 = M3
                  IF M3 = 0 THEN
                   CALL threej(R0, AJ1, 6!, AJ1, -AM1, AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1)) * R0 * RMC(3) * B6Q(M3)
                   ELSE
                   CALL threej(R0, AJ1, 6!, AJ1, -AM1, AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1)) * R0 * RMC(3) * B6Q(M3)
                   CALL threej(R0, AJ1, 6!, AJ1, -AM1, -AM3, AM2)
                   EMAT(IA1, IA2) = EMAT(IA1, IA2) + ((-1) ^ (AJ1 - AM1 + M3)) * R0 * RMC(3) * B6Q(M3)
                  END IF
              NEXT M3
          NEXT AM2
      NEXT AM1

REM   FIND THE EIGENVALUES OF THE MATRIX EMAT
REM   ***
REM    * FIRST DEFINE THE MATRIX
      N = 2 * AJ1 + 1
      FOR I = 1 TO N
          FOR J = 1 TO N
              A(I, J) = EMAT(I, J)
          NEXT J
      NEXT I
REM    * EIGENVALUE ANALYSIS BEGINS
      PRINT "FINDING THE EIGENVALUES"
      FOR IP = 1 TO N
        FOR IQ = 1 TO N
          V(IP, IQ) = 0!
        NEXT IQ
        V(IP, IP) = 1!
      NEXT IP
      FOR IP = 1 TO N
        B(IP) = A(IP, IP)
        D(IP) = B(IP)
        Z(IP) = 0!
      NEXT IP
      NROT = 0
      FOR I = 1 TO 50
        SM = 0!
        FOR IP = 1 TO N - 1
          FOR IQ = IP + 1 TO N
            SM = SM + ABS(A(IP, IQ))
          NEXT IQ
        NEXT IP
        IF SM = 0! THEN GOTO RET:
        IF I < 4 THEN
          TRESH = .2 * SM / N ^ 2
        ELSE
          TRESH = 0!
        END IF
        FOR IP = 1 TO N - 1
          FOR IQ = IP + 1 TO N
            G = 100! * ABS(A(IP, IQ))
     IF (I > 4) AND (ABS(D(IP)) + G = ABS(D(IP))) AND (ABS(D(IQ)) + G = ABS(D(IQ))) THEN
                  A(IP, IQ) = 0!
     ELSE
          IF ABS(A(IP, IQ)) > TRESH THEN
              H = D(IQ) - D(IP)
            IF (ABS(H) + G = ABS(H)) THEN
                T = A(IP, IQ) / H
            ELSE
                THETA = .5 * H / A(IP, IQ)
                T = 1! / (ABS(THETA) + SQR(1! + THETA ^ 2))
                IF THETA < 0! THEN
                T = -T
                END IF
             END IF
              C = 1! / SQR(1 + T ^ 2)
              S = T * C
              TAU = S / (1! + C)
              H = T * A(IP, IQ)
              Z(IP) = Z(IP) - H
              Z(IQ) = Z(IQ) + H
              D(IP) = D(IP) - H
              D(IQ) = D(IQ) + H
              A(IP, IQ) = 0!
              FOR J = 1 TO IP - 1
                G = A(J, IP)
                H = A(J, IQ)
                A(J, IP) = G - S * (H + G * TAU)
                A(J, IQ) = H + S * (G - H * TAU)
              NEXT J
              FOR J = IP + 1 TO IQ - 1
                G = A(IP, J)
                H = A(J, IQ)
                A(IP, J) = G - S * (H + G * TAU)
                A(J, IQ) = H + S * (G - H * TAU)
              NEXT J
              FOR J = IQ + 1 TO N
                G = A(IP, J)
                H = A(IQ, J)
                A(IP, J) = G - S * (H + G * TAU)
                A(IQ, J) = H + S * (G - H * TAU)
              NEXT J
              FOR J = 1 TO N
                G = V(J, IP)
                H = V(J, IQ)
                V(J, IP) = G - S * (H + G * TAU)
                V(J, IQ) = H + S * (G - H * TAU)
              NEXT J
              NROT = NROT + 1
            END IF
            END IF
        NEXT IQ
      NEXT IP
        FOR IP = 1 TO N
          B(IP) = B(IP) + Z(IP)
          D(IP) = B(IP)
          Z(IP) = 0!
      NEXT IP
    NEXT I
     PRINT "50 iterations should never happen"
     PRINT "PROGRAM ABORTED IN EIGENVALUE ANALYSIS"
RET:

REM **
REM   NAME OF OUTPUT FILES OF ENERGY LEVELS
      PRINT "EIGENVALUES AND EIGENVECTORS HAVE BEEN DETERMINED"
      PRINT "  "
      PRINT "OUTPUT WILL BE IN A SINGLE FILE CONTAINING"
      PRINT "EIGENVALUES IN ARBITRARY ORDER.  DO YOU WISH"
      PRINT "TO OUTPUT A SECOND FILE CONTAINING BOTH "
      INPUT "EIGENVECTORS AND EIGENVALUES (Y/y FOR YES) ="; ANSO$

      IF ANSO$ = "y" OR ANSO$ = "Y" THEN
       INPUT "NAME OF EIGENVALUE FILE = "; FILE2$
       INPUT "NAME OF EIGENVECTOR AND EIGENVALUE FILE = "; FILE3$
       PRINT "  "
       PRINT "OUTPUT FILE FOR EIGENVALUES = "; FILE2$
       PRINT "OUTPUT FILE FOR EIGENVECTORS AND EIGENVALUES = "; FILE3$
       ELSE
       INPUT "NAME OF EIGENVALUE FILE = "; FILE2$
       PRINT "  "
       PRINT "OUTPUT FILE FOR EIGENVALUES = "; FILE2$
      END IF
     
        OPEN FILE2$ FOR OUTPUT AS #8
      IF ANSO$ = "Y" OR ANSO$ = "y" THEN
        OPEN FILE3$ FOR OUTPUT AS #9
      END IF
       IK = 2 * AJ1 + 1
       FOR I = 1 TO IK
        IF ANSO$ = "Y" OR ANSO$ = "y" THEN
           PRINT #9, "ENERGY = "; D(I)
           PRINT #9, "WITH EIGENVECTOR:"
           FOR J = 1 TO IK
             PRINT #9, V(J, I)
           NEXT J
             PRINT #8, USING "########.#"; D(I)
           ELSE
             PRINT #8, USING "########.#"; D(I)
           END IF
       NEXT I

       PRINT "PROGRAM RUN IS COMPLETED SUCCESSFULLY"
    
       END

DEFSNG I-J, M-N
SUB threej (R0, AJ1, AJ2, AJ3, AM1, AM2, AM3)
       DEFINT J, M-N
       DEFDBL F, R
       DIM F(0 TO 30)

REM DEFINE THE FACTORIAL FUNCTION
       FOR I = 0 TO 30
           F(I) = 1
           FOR J = 1 TO I
             F(I) = F(I) * J
           NEXT J
       NEXT I

REM   SELECTION RULES
     IF AM1 + AM2 + AM3 <> 0 THEN GOTO ZERO:
     IF AJ1 + AJ2 - AJ3 < 0 THEN GOTO ZERO:
     IF AJ3 + AJ1 - AJ2 < 0 THEN GOTO ZERO:
     IF AJ3 + AJ2 - AJ1 < 0 THEN GOTO ZERO:
     IF AJ1 + AJ2 + AJ3 + 1 < 0 THEN GOTO ZERO:
     IF AJ1 - AM1 < 0 THEN GOTO ZERO:
     IF AJ2 - AM2 < 0 THEN GOTO ZERO:
     IF AJ3 - AM3 < 0 THEN GOTO ZERO:
     IF AJ1 + AM1 < 0 THEN GOTO ZERO:
     IF AJ2 + AM2 < 0 THEN GOTO ZERO:
     IF AJ3 + AM3 < 0 THEN GOTO ZERO:

   R4 = F(CINT(AJ1 + AJ2 - AJ3)) * F(CINT(AJ1 - AM1)) * F(CINT(AJ2 - AM2)) * F(CINT(AJ3 - AM3)) * F(CINT(AJ3 + AM3))
   R5 = F(CINT(AJ1 + AJ2 + AJ3 + 1)) * F(CINT(AJ3 + AJ1 - AJ2)) * F(CINT(AJ3 + AJ2 - AJ1)) * F(CINT(AJ1 + AM1)) * F(CINT(AJ2 + AM2))
   R6 = 0
   FOR J7 = 0 TO 25
       IF AJ1 + AM1 + J7 < 0 THEN GOTO J7R:
       IF AJ2 + AJ3 - AM1 - J7 < 0 THEN GOTO J7R:
       IF AJ3 + AM3 - J7 < 0 THEN GOTO J7R:
       IF AJ1 - AM1 - J7 < 0 THEN GOTO J7R:
       IF AJ2 - AJ3 + AM1 + J7 < 0 THEN GOTO J7R:
       R8 = F(CINT(AJ1 + AM1 + J7)) * F(CINT(AJ2 + AJ3 - AM1 - J7)) * (-1) ^ (CINT(AJ1 - AM1 - J7))
       R9 = F(J7) * F(CINT(AJ3 + AM3 - J7)) * F(CINT(AJ1 - AM1 - J7)) * F(CINT(AJ2 - AJ3 + AM1 + J7))
       R6 = R6 + R8 / R9
J7R:  NEXT J7
       R0 = SQR(R4 / R5) * R6 * (-1) ^ CINT((AJ1 - AJ2 - AM3))
       GOTO FIN:
ZERO:  R0 = 0
FIN: END SUB