PDP-10 Archive: decus/20-0077/nmrsim.for from decus_20tap3

Trailing-Edge - PDP-10 Archives - decus_20tap3_198111 - decus/20-0077/nmrsim.for

There is 1 other file named nmrsim.for in the archive. Click here to see a list.

C	NMR SIMULATION AND PLOTTING PROGRAM
C
C
C
C	9 JUNE 1973
C	JAMES S. EVANS
C
C
C	PROGRAM IS BASED ON FORMALISM, THEOREMS, AND METHODS
C	DETAILED BY P. L. CORIO, "STRUCTURE OF HIGH-RESOLUTION
C	NMR SPECTRA" (ACADEMIC PRESS, 1966)
C
	IMPLICIT INTEGER (P,W)
C
C	AS	SYMBOLS FOR AVAILABLE NUCLEI
C	AX	NUCLEAR SPIN QUANTUM NUMBERS
C	AG	GYROMAGNETIC RATIOS (RECIPROCAL MILLIGAUSS-SEC)
C
	COMMON /NUCL/ AS(10),AX(10),AG(10)
C
C	A	HAMILTONIAN MATRIX
C	S	MATRIX OF EIGENVECTORS (STORED COLUMNWISE)
C	T	TRANSITION PROBABILITY MATRIX
C	R	TEMPORARY STORAGE FOR HAMIL
C	X	TEMPORARY STORAGE FOR SPINS AND EIGEN2
C	MY	TEMPORARY STORAGE FOR EIGEN2
C	P	POINTING VECTOR (AVOIDS SORTING SPIN FUNCTIONS)
C	Z	SPIN FUNCTIONS (Z-COMPONENT FOR EACH NUCLEUS,
C		TOTAL Z-COMPONENT IN COL 0, ONE ROW PER SPIN STATE)
C	GZ	GYROMAGNETIC RATIOS OF NUCLEI IN MOLECULE
C	YJ	CHEMICAL SHIFTS (ON-DIAGONAL) AND COUPLING CONSTANTS
C		(OFF-DIAGONAL)
C	YL	SPINS OF NUCLEI IN MOLECULE
C	BS	SYMBOLS OF NUCLEI IN MOLECULE
C	W	SIZE OF HAMILTONIAN MATRIX FOR EACH SPIN
C	WV	TRANSFER VECTOR (USED WITH P TO ACCESS SPIN FUNCTIONS)
C
	COMMON /CALC/ A(35,35),S(35,35),T(35,35),R(35),X(35),MY(35),
	1 P(128),Z(128,0/7),GZ(7),YJ(7,7),YL(7),BS(0/7),W(15),WV(15)
C
C	E	FREQUENCIES IN LINE SPECTRUM
C	YI	INTENSITIES IN LINE SPECTRUM (REAL I IN SUBS.
C		ORDER AND SHAPE)
C	IX	FREQUENCY SCALE FOR PLOT (INTEGER X IN SUBS.
C		ORDER AND SHAPE)
C	Y	INTENSITY SCALE FOR PLOT
C
	COMMON /SPEC/ E(500),YI(500),IX(500),Y(0/510)
C
C	IO2,IO3,IO4 INPUT,OUTPUT TO DISK (OR MAG TAPE)
C	IO5,IO6 INPUT,OUTPUT TO USER'S TERMINAL
C
	COMMON /IODEV/ IO2,IO3,IO4,IO5,IO6
C
C	KWIT	SIGNAL FOR ESCAPE KEY STRUCK
C	L	PEN UP (.NE. 0) OR DOWN (.EQ. 0)
C	X0	MINIMUM X BOUNDARY
C	X1	MAXIMUM X BOUNDARY
C	X2	NEW X
C	X3	OLD X, ALSO TOTAL X DISPLACEMENT
C	X4	INTERMEDIATE X DISPLACEMENT
C	X5	DIGITAL RESOLUTION OF TSP-12 PLOTTER INTERFACE
C	Y0-Y5	DEFINED SIMILARLY
C
	COMMON /PLOT1/ KWIT,L,X0,X1,X2,X3,X4,X5,Y0,Y1,Y2,Y3,Y4,Y5
C
	INTEGER X2,X3,X4,X5,Y2,Y3,Y4,Y5
C
	DATA AS/'H','D','B10','B11','C13','N14','N15','F','P','LAST'/
	DATA AX/0.5, 1.0, 3.0, 1.5, 0.5, 1.0, 0.5, 0.5, 0.5, 0.0/
	DATA AG/26.7512, 4.10641, 2.87465, 8.58276, 6.72597,
	1 1.93289, 2.71123, 25.1668, 10.829, 0./
C
C	PROGRAM INITIALIZATION
C
	IO2=20
	IO3=21
	IO4=22
	IO5=5
	IO6=5
	BS(0)=' '
	N9=500
	R8=1.0E-5
	X5=510
	Y5=510
C
	WRITE(IO6,6000)
	NCHK=0
	GO TO 20

C
C	SUPERVISORY SECTION
C
C	MODE	CONTROLS FLOW OF OPERATIONS
C	NCHK	PREVENTS ACCESS TO ROUTINES LACKING NEEDED DATA
C	NGO	SPECIFIES DESTINATION, ALSO ERROR RETURN
C
C	ERROR REENTRY AT STATEMENT 10
C
10	WRITE(IO6,6001) NGO
C
C	NORMAL REENTRY AT STATEMENT 20
C
20	WRITE(IO6,6002)
	READ(IO5,5000) NGO
	IF(NGO .GE. 0 .AND. NGO .LE. 30) GO TO 50
C
C	NONEXISTENT COMMAND CODES TRAP TO STATEMENT 30
C
30	WRITE(IO6,6003)
	GO TO 20
C
C
C	ROUTING TO VARIOUS OPTIONS
C
50	KWIT=0
	GO TO (90,100,200,30,30,30,30,30,30,30,
	1 1000,1100,1200,1300,1400,1500,1600,30,30,30,
	2 30,2100,2200,2300,2400,2500,2600,2700,2800,30,3000),NGO+1
C
C	HELP MESSAGE
C
90	WRITE(IO6,6004)
C
C	CHANGE TO SEQUENTIAL RUN MODE (FOR DEBUGGING)
C
100	MODE=1
	GO TO 20
C
C	CHANGE TO EPISODIC RUN MODE (FOR DEBUGGING)
C
200	MODE=2
	GO TO 20

C
C	CALCULATION OF LINE SPECTRUM (OPTIONS 10-16)
C
C
C	LIST AVAILABLE NUCLEI (OPTION 10)
C
1000	WRITE(IO6,6005)
	I=1
1010	IF(AS(I) .EQ. 'LAST') GO TO 1020
	WRITE(IO6,6006) AS(I)
	I=I+1
	GO TO 1010
1020	WRITE(IO6,6007)
	NCHK=NCHK .OR. "2000
	GO TO 20
C
C	DEFINE NEW MOLECULE (OPTION 11)
C
C	YI5	FACTOR CONTRIBUTING TO THEORETICAL INTENSITY
C	N8	NUMBER OF NUCLEI IN MOLECULE WITH NONZERO SPIN
C	YL9	SUM OF ALL SPINS IN MOLECULE (HIGHEST SPIN STATE)
C
1100	YI5=2./3.
	YL9=0.
	MODE=1
1110	NCHK=NCHK .AND. "777774003777
	WRITE(IO6,6008)
	READ(IO5,5000) N8
	IF(N8 .LE. 0 .OR. N8 .GT. 7) GO TO 1110
	WRITE(IO6,6009)
	DO 1150 I=1,N8
1120	WRITE(IO6,6010) I
	READ(IO5,5001) BS(I)
	J=0
1130	J=J+1
	IF(AS(J) .NE. 'LAST') GO TO 1140
	WRITE(IO6,6011) BS(I)
	GO TO 1120
1140	IF(AS(J) .NE. BS(I)) GO TO 1130
	GZ(I)=AG(J)
	YL(I)=AX(J)
	X(I)=AX(J)
	YL9=YL9+AX(J)
	YI5=YI5*(2.*AX(J)+1.)
1150	CONTINUE
	CALL SPINS (YL9,N8)
	IF(YL9 .GT. 0.) GO TO 1160
	WRITE(IO6,6012)
	GO TO 1100
1160	NCHK=NCHK .OR. "4000
C
C	ENTER COUPLING CONSTANTS (OPTION 12)
C
1200	IF((NCHK .AND. "4000) .EQ. "4000) GO TO 1210
	NGO=12
	GO TO 10
1210	NCHK=NCHK .AND. "777777567777
	WRITE(IO6,6013)
	DO 1220 I=1,N8-1
	DO 1220 J=I+1,N8
	WRITE(IO6,6014) I,J
	READ(IO5,5002) YJ(I,J)
	YJ(J,I)=YJ(I,J)
1220	CONTINUE
	NCHK=NCHK .OR. "10000
	IF(MODE .EQ. 2) GO TO 20
C
C	ENTER CHEMICAL SHIFTS (OPTION 13)
C
1300	IF((NCHK .AND. "4000) .EQ. "4000) GO TO 1310
	NGO=13
	GO TO 10
1310	NCHK=NCHK .AND. "777777557777
	WRITE(IO6,6015)
	DO 1320 I=1,N8
	IF(BS(I) .NE. BS(I-1)) WRITE(IO6,6007)
	WRITE(IO6,6016) I,BS(I)
	READ(IO5,5002) YJ(I,I)
1320	CONTINUE
	NCHK=NCHK .OR. "20000
	IF(MODE .EQ. 2) GO TO 20
C
C	ENTER NUCLEUS OBSERVED (OPTION 14)
C
C	ASO	SYMBOL FOR NUCLEUS IN RESONANCE
C	YI9	TOTAL THEORETICAL INTENSITY
C		(CORIO, EQS. 2.12 AND 2.13, P. 288)
C
1400	IF((NCHK .AND. "4000) .EQ. "4000) GO TO 1410
	NGO=14
	GO TO 10
1410	NCHK=NCHK .AND. "777777537777
	WRITE(IO6,6017)
	READ(IO5,5001) ASO
	YI9=0.0
	DO 1420 I=1,N8
	IF(BS(I) .NE. ASO) GO TO 1420
	G=GZ(I)
	YI9=YI9+YL(I)*(YL(I)+1.)
1420	CONTINUE
	IF(YI9 .GT. 0.0) GO TO 1430
	WRITE(IO6,6018) ASO
	NGO=14
	GO TO 10
1430	YI9=YI5*YI9
	NCHK=NCHK .OR. "40000
	IF(MODE .EQ. 2) GO TO 20
C
C	ENTER INTENSITY THRESHOLD (OPTION 15)
C
C	YI2	CUT-OFF INTENSITY TO REJECT NUMERICAL ARTIFACTS
C
1500	IF((NCHK .AND. "44000) .EQ. "44000) GO TO 1510
	NGO=15
	GO TO 10
1510	NCHK=NCHK .AND. "777777477777
	WRITE(IO6,6019) YI9
	READ(IO5,5002) YI2
	YI2=SQRT(ABS(YI2))
	NCHK=NCHK .OR. "100000
	IF(MODE .EQ. 2) GO TO 20
C
C	START DIAGONALIZATIONS (OPTION 16)
C
C	FNAME	FILENAME ON THE DISK WILL BE FNAME.DAT
C	YI4	TOTAL CALCULATED INTENSITY
C	N7	TOTAL NUMBER OF UNSORTED AND UNCOMBINED LINES
C	E1	MOST NEGATIVE FREQUENCY IN LINE SPECTRUM
C	E2	MOST POSITIVE FREQUENCY IN LINE SPECTRUM
C
1600	IF((NCHK .AND. "174000) .EQ. "174000) GO TO 1610
	NGO=16
	GO TO 10
1610	NCHK=NCHK .AND. "777777577777
	WRITE(IO6,6020)
	READ(IO5,5001) FNAME
	CALL OFILE(IO3,FNAME)
	WRITE(IO3,6021) G
	CALL HAMIL (ASO,R8,YL9,N8,YI2,YI4,N7,E1,E2,KWIT)
	IF(KWIT .NE. 32) GO TO 1620
	WRITE(IO6,6022) FNAME
	GO TO 1630
1620	END FILE IO3
	WRITE(IO6,6023) FNAME,YI4,N7,E2,E1
1630	REWIND IO3
	IF(MODE .EQ. 1) MODE=2
	NCHK=NCHK .OR. "200000
	GO TO 20

C	PLOTTING SECTION (OPTIONS 21-28)
C
C
C	RETRIEVE LINE SPECTRUM FOR PLOTTING (OPTION 21)
C
C	FILE	FILENAME TO BE PLOTTED (WON'T INTERFERE WITH FNAME)
C	GP	GYROMAGNETIC RATIO
C
2100	NCHK=NCHK .AND. "7777777
	WRITE(IO6,6020)
	READ(IO5,5001) FILE
	CALL IFILE(IO2,FILE)
	READ(IO2,6021) GP
	NCHK=NCHK .OR. "10000000
C
C	ENTER MINIMUM RESOLUTION (OPTION 22)
C
C	R1	RESOLUTION (HZ)
C	R0	RESOLUTION (SEC)
C	N7	NUMBER OF SORTED AND COMBINED LINES
C	N9	DIMENSIONED SIZE OF ENERGY AND INTENSITY ARRAYS
C
2200	IF(MODE .EQ. 2) MODE=1
	IF((NCHK .AND. "10000000) .EQ. "10000000) GO TO 2210
	NGO=22
	GO TO 10
2210	NCHK=NCHK .AND. "770017777777
	WRITE(IO6,6024)
	READ(IO5,5002) R1
	R0=1./(R1+1.0E-10)
	CALL ORDER (N7,N9,R0,KWIT)
	IF(KWIT .GT. 0) GO TO 2220
	R1=1./R0
	WRITE(IO6,6025) N9
	WRITE(IO6,6026) R1
2220	WRITE(IO6,6027) N7
	WRITE(IO6,6031) E(N7),E(1)
	NCHK=NCHK .OR. "20000000
	IF(MODE .EQ. 2) GO TO 20
C
C	TABULATE LINE SPECTRUM (OPTION 23)
C
2300	IF((NCHK .AND. "30000000) .EQ. "30000000) GO TO 2310
	NGO=23
	GO TO 10
2310	NCHK=NCHK .AND. "777737777777
	WRITE(IO6,6028)
	READ(IO5,5002) J
	IF((J .NE. 1) .AND. (J .NE. 2)) GO TO 2360
	WRITE(IO6,6038)
	READ(IO5,5001) FILE2
	CALL OFILE(IO4,FILE2)
	WRITE(IO4,6007)
	WRITE(IO4,5001) FILE2
	WRITE(IO4,6007)
	WRITE(IO4,6029)
	WRITE(IO4,6030)(E(I),YI(I),I=1,N7)
	WRITE(IO4,6007)
	WRITE(IO4,6027) N7
	WRITE(IO4,6026) R1
	END FILE IO4
	REWIND IO4
	IF(J .NE. 2) GO TO 2350
	WRITE(IO6,6039)
	READ(IO5,5002) YI3
	IF(YI3 .EQ. 0.) YI3=1.0E-10
	WRITE(IO6,6029)
	KWIT=0
	DO 2320 I=1,N7
	IF(YI(I) .GE. YI3) WRITE(IO6,6030) E(I),YI(I)
	CALL TTYTST(KWIT,K)
	IF(KWIT .EQ. 32) GO TO 2360
2320	CONTINUE
2350	WRITE(IO6,6007)
	WRITE(IO6,6027) N7
	WRITE(IO6,6026) R1
	WRITE(IO6,6007)
2360	NCHK=NCHK .OR. "40000000
	IF(MODE .EQ. 2) GO TO 20
C
C	ENTER X-SCALE (OPTION 24)
C
2400	IF((NCHK .AND. "30000000) .EQ. "30000000) GO TO 2410
	NGO=24
	GO TO 10
2410	NCHK=NCHK .AND. "771277777777
	WRITE(IO6,6031) E(N7),E(1)
	WRITE(IO6,6032)
	READ(IO5,5002) X0,X1
	NCHK=NCHK .OR. "100000000
	IF(MODE .EQ. 2) GO TO 20
C
C	ENTER H1,T1,T2
C
C	H1	RF POWER (MILLIGAUSS)
C	T1,T2	RELAXATION TIMES (SEC)
C
2500	IF((NCHK .AND. "30000000) .EQ. "30000000) GO TO 2510
	NGO=25
	GO TO 10
2510	NCHK=NCHK .AND. "771177777777
	WRITE(IO6,6033)
	READ(IO5,5002) H1
	WRITE(IO6,6034)
	READ(IO5,5002) T1,T2
	NCHK=NCHK .OR. "200000000
	IF(MODE .EQ. 2) GO TO 20
C
C	ENTER Y-SCALE (OPTION 26)
C
C	Y9	HEIGHT OF BIGGEST PEAK TO BE PLOTTED
C
2600	IF((NCHK .AND. "330000000) .EQ. "330000000) GO TO 2610
	NGO=26
	GO TO 10
2610	NCHK=NCHK .AND. "771377777777
	CALL SHAPE (N7,GP,H1,T1,T2,Y9)
	WRITE(IO6,6035) Y9
	READ(IO5,5002) Y0,Y1
	NCHK=NCHK .OR. "400000000
	IF(MODE .EQ. 2) GO TO 20
C
C	ENTER PLOT SPEED (OPTION 27)
C
C	P4	NUMBER OF BITS (OUT OF 511) FOR TYPICAL PEN MOVEMENT
C
2700	NCHK=NCHK .AND. "776777777777
	WRITE(IO6,6036) X5
	READ(IO5,5000) P4
	IF(P4 .LT. 1) P4=1
	NCHK=NCHK .OR. "1000000000
	IF(MODE .EQ. 2) GO TO 20
C
C	PLOT WITH CURRENT PARAMETERS (OPTION 28)
C
2800	IF((NCHK .AND. "1730000000) .EQ. "1730000000) GO TO 2810
	NGO=28
	GO TO 10
2810	NCHK=NCHK .AND. "775777777777
	CALL PCTRL(1)
	L=1
	J=0
	YY=-Y(0)-1
	DO 2840 X2=0,X5
	IF(Y(X2) .NE. YY) GO TO 2830
	J=J+1
	IF(J .LT. P4) GO TO 2840
2830	J=0
	YY=Y(X2)
	CALL YPLOT (YY,P4)
	IF(KWIT .EQ. 32) GO TO 2850
2840	CONTINUE
	CALL PCTRL(3)
2850	IF(MODE .EQ. 1) MODE=2
	NCHK=NCHK .OR. "2000000000
	GO TO 20
C
C	EXIT FROM PROGRAM TO SYSTEM MONITOR
C
3000	CONTINUE
	WRITE(IO6,6037)
	CALL EXIT

C
5000	FORMAT(I)
5001	FORMAT(A5)
5002	FORMAT(2G)
C
6000	FORMAT(/' NMR SIMULATION AND PLOTTING PROGRAM'/
	1 '0USER MUST GIVE TWO-DIGIT REPLY WHENEVER PROGRAM'/
	1 ' ASKS "NEXT?"  TO EXIT, TYPE 30.  FOR HELP, TYPE 0.')
6001	FORMAT(' COMMAND SEQUENCE ERROR AT OPTION'I3,'--BACK UP?')
6002	FORMAT('0NEXT? '$)
6003	FORMAT(' COMMAND CODE NOT RECOGNIZED')
6004	FORMAT('  0=TYPE THIS TEXT'//
	1 ' 10=LIST AVAILABLE NUCLEI'/
	1 ' 11=DEFINE NEW SPIN SYSTEM'/
	1 ' 12=ENTER COUPLING CONSTANTS'/
	1 ' 13=ENTER CHEMICAL SHIFTS'/
	1 ' 14=ENTER NUCLEUS OBSERVED'/
	1 ' 15=ENTER INTENSITY THRESHOLD'/
	1 ' 16=START DIAGONALIZATIONS'//
	1 ' 21=RETRIEVE LINE SPECTRUM TO PLOT'/
	1 ' 22=ENTER MINIMUM RESOLUTION'/
	1 ' 23=TABULATE LINE SPECTRUM'/
	1 ' 24=ENTER X-SCALE'/
	1 ' 25=ENTER H1,T1,T2'/
	1 ' 26=ENTER Y-SCALE'/
	1 ' 27=ENTER PLOT SPEED'/
	1 ' 28=PLOT WITH CURRENT PARAMETERS'//
	1 ' 30=QUIT')
6005	FORMAT(' AVAILABLE NUCLEI: '$)
6006	FORMAT('+'A5,$)
6007	FORMAT(' ')
6008	FORMAT(' HOW MANY NUCLEI (MAX=7)? '$)
6009	FORMAT(' IDENTIFY NUCLEI (ALL OF ONE TYPE,'
	1 ' THEN NEXT TYPE, ETC.)'/)
6010	FORMAT('+#'I1,'  '$)
6011	FORMAT(' 'A5,' NOT AVAILABLE: REENTER'$)
6012	FORMAT(' TOO MANY SPIN FUNCTIONS: USE FEWER NUCLEI')
6013	FORMAT(' PAIR    COUPLING CONSTANT (HZ)'/)
6014	FORMAT('+'2I2,'      '$)
6015	FORMAT(' ATOM#  TYPE  CHEMICAL SHIFT (HZ)')
6016	FORMAT('+'I3,5X,A5,1X,$)
6017	FORMAT(' TYPE OF NUCLEUS OBSERVED? '$)
6018	FORMAT(' NUCLEUS 'A5,' NOT DEFINED IN MOLECULE')
6019	FORMAT(' TOTAL INTENSITY EXPECTED = 'F6.1/
	1 ' MINIMUM INTENSITY? '$)
6020	FORMAT(' NAME FOR LINE SPECTRUM DISK FILE (MAX=5 CHARS)? '$)
6021	FORMAT(E15.8)
6022	FORMAT(' DIAGONALIZATION ABORTED, CAN REUSE FILE 'A5)
6023	FORMAT(' DIAGONALIZATIONS COMPLETED FOR FILE 'A5/
	1 ' TOTAL COMPUTED INTENSITY = 'F10.5,' IN 'I5,' LINES'/
	1 ' FROM'F10.3,' TO 'F10.3,' HZ')
6024	FORMAT(' MINIMUM MEANINGFUL RESOLUTION (HZ)? '$)
6025	FORMAT(' MORE THAN 'I4,' LINES WITH DEGRADED')
6026	FORMAT(' RESOLUTION = 'F8.4,' HZ')
6027	FORMAT(' NUMBER OF LINES = 'I4)
6028	FORMAT(' FOR LISTING OF LINES, SPECIFY 1=DSK, 2=TTY ALSO,',
	1 ' ELSE SKIP? '$)
6029	FORMAT('      FREQUENCY      INTENSITY'/)
6030	FORMAT(2F15.5)
6031	FORMAT(' LINES EXTEND FROM'F10.3,' TO'F10.3,' HZ')
6032	FORMAT(' WHAT LEFT,RIGHT FREQUENCIES (HZ) FOR X-AXIS? '$)
6033	FORMAT(' WHAT RF POWER (MILLIGAUSS,APPROX 0.01)? '$)
6034	FORMAT(' WHAT T1,T2 (SEC)? '$)
6035	FORMAT(' MAXIMUM INTENSITY (ARBITRARY) ='F8.4/
	1 ' WHAT BOTTOM,TOP INTENSITIES FOR Y AXIS? '$)
6036	FORMAT(' WHAT PLOT SPEED (PARTS IN'I4,')? '$)
6037	FORMAT(' DISK FILES CREATED BY THIS PROGRAM HAVE .DAT'/
	1 ' EXTENSION AND 057 PROTECTION CODE')
6038	FORMAT(' NAME FOR LISTING FILE (MAX=5 CHARS)? '$)
6039	FORMAT(' FOR TTY, SPECIFY 0=TYPE ALL LINES,',
	1 ' ELSE ENTER MIN INTENSITY? '$)
	END

	SUBROUTINE SPINS (YL9,N8)
C
	IMPLICIT INTEGER (P,W)
C
C	PREPARES SPIN BASIS FUNCTIONS FOR THE MOLECULE
C
	COMMON /CALC/ A(35,35),S(35,35),T(35,35),R(35),X(35),MY(35),
	1 P(128),Z(128,0/7),GZ(7),YJ(7,7),YL(7),BS(0/7),W(15),WV(15)
C
C
C	W9	TOTAL NUMBER OF HAMILTONIAN MATRICES
C	WZ9	TOTAL NUMBER OF SPIN MICROSTATES
C
	W9=2.*YL9+1.
	DO 110 I=1,W9
110	W(I)=0
C
	DO 140 I=1,128
	Z(I,0)=0.
	DO 120 J=1,N8
	Z(I,J)=X(J)
120	Z(I,0)=Z(I,0)+X(J)
	K=IFIX(YL9-Z(I,0))+1
	W(K)=W(K)+1
	L=1
130	X(L)=X(L)-1.
	IF(YL(L) .GE. -X(L)) GO TO 140
	X(L)=YL(L)
	L=L+1
	IF(L .GT. N8) GO TO 150
	GO TO 130
140	CONTINUE
	YL9=-YL9
	RETURN
C
150	WZ9=I
	WV(1)=1
	DO 160 I=2,W9
160	WV(I)=WV(I-1)+W(I)
C
	DO 170 I=1,WZ9
	L=IFIX(YL9-Z(I,0))+1
	P(WV(L))=I
170	WV(L)=WV(L)-1
C
	WV(1)=1
	DO 180 I=2,W9
180	WV(I)=WV(I-1)+W(I)
	RETURN
	END

	SUBROUTINE HAMIL (ASO,R8,L9,N8,I2,I4,N7,E1,E2,ESCAP)
C
C
C	PREPARES HAMILTONIAN MATRICES, DIAGONALIZES THEM,
C	COMPUTES TRANSITION ENERGIES AND INTENSITIES
C
	COMMON /CALC/ A(35,35),S(35,35),T(35,35),R(35),X(35),MY(35),
	1 P(128),Z(128,0/7),GZ(7),YJ(7,7),YL(7),BS(0/7),W(15),WV(15)
C
	COMMON /IODEV/ IO2,IO3,IO4,IO5,IO6
C
	INTEGER ESCAP,P,W,WV
	REAL I2,I4,I7,I8,L9
C
C
C	H	SPIN VALUE FOR CURRENT MATRIX
C	L9	MAX SPIN (YL9 IN MAIN PROGRAM)
C	N	SIZE OF MATRIX
C	I2	INTENSITY CUT-OFF (YI2 IN MAIN PROGRAM)
C	I4	TOTAL INTENSITY (YI4 IN MAIN PROGRAM)
C	R8	CONVERGENCE CRITERION FOR DIAGONALIZATION
C	I7	INTENSITY OF LINE
C	E7	ENERGY OF LINE
C	N7	COUNTER FOR NUMBER OF LINES
C
	I4=0.
	N7=0
	E1=1.0E10
	E2=-E1
	H=L9
C
100	IH=L9-H+1.
	N=W(IH)
	I0=WV(IH)-N
	IF(H .EQ. L9) GO TO 200
C
C	MAKE ALLOWED TRANSITION MATRIX T FROM OLD EIGENVECTORS S
C	(THIS SECTION SKIPPED FOR H=L9)
C
C	MATRIX T IS CONSTRUCTED SO THAT ELEMENTS OF THE MATRIX
C	PRODUCT T.S WILL GIVE INTENSITIES (CORIO, EQ. 3.41, P. 166)
C
C	SELECTION RULES IN THE X APPROXIMATION (CORIO, PP. 286-7)
C		(1) SPIN DECREASE OF UNITY FOR ASO
C		(2) OTHER NUCLEI MAY NOT CHANGE
C		(3) ONLY A SINGLE ASO-TYPE NUCLEUS MAY CHANGE
C
	DO 110 I=1,N6
	R(I)=A(I,I)
	DO 110 J=1,N
	T(I,J)=0.
110	CONTINUE
C
	DO 170 I=1,N
	I1=P(I0+I)
	DO 160 J=1,N6
	J1=P(I6+J)
	L=0
	DO 140 K=1,N8
	I8=Z(I1,K)-Z(J1,K)
	IF(I8 .EQ. 0.) GO TO 140
	IF((I8 .NE. -1.) .OR. (BS(K) .NE. ASO)) GO TO 160
	K1=K
	L=L+1
140	CONTINUE
	IF(L .NE. 1) GO TO 160
	C=SQRT((YL(K1)+Z(J1,K1))*(YL(K1)-Z(J1,K1)+1.))
	DO 150 L=1,N6
150	T(L,I)=T(L,I)+C*S(J,L)
160	CONTINUE
170	CONTINUE
C
C	MAKE HAMILTONIAN SUBMATRIX FOR SPIN H
C
C	SIGNS OF MATRIX ELEMENTS ARE OPPOSITE TO THOSE OF CORIO,
C	BUT COMPENSATION OCCURS WHEN TRANSITION ENERGY IS COMPUTED
C	BY DIFFERENCE
C
200	DO 280 I=1,N
	I1=P(I0+I)
	A(I,I)=0.
C
C	ON-DIAGONAL ELEMENTS (CORIO, EQ. 2.15, P. 152)
C
	DO 220 J=1,N8
	A(I,I)=A(I,I)+YJ(J,J)*Z(I1,J)
	IF(J .EQ. N8) GO TO 220
	DO 210 K=J+1,N8
210	A(I,I)=A(I,I)+YJ(J,K)*Z(I1,J)*Z(I1,K)
220	CONTINUE
C
C	OFF-DIAGONAL ELEMENTS (CORIO, EQ. 2.16, P. 153)
C
C		I8=SPIN CHANGE FOR NUCLEUS K
C		SPIN-SPIN INTERACTION INCREASES SPIN OF NUCLEUS J1
C		 AND DECREASES SPIN OF K1 BY ONE UNIT EACH
C		NO OTHER NUCLEUS MAY CHANGE
C
	IF(I .EQ. N) GO TO 280
	DO 270 J=I+1,N
	L=0
	A(I,J)=0.
	A(J,I)=0.
C
	DO 260 K=1,N8
	I8=Z(I1,K)-Z(P(I0+J),K)
	IF(ABS(I8) .GT. 1.) GO TO 270
	IF(I8) 230,260,240
230	J1=K
	GO TO 250
240	K1=K
250	L=L+1
	IF(L .GT. 2) GO TO 270
260	CONTINUE
C
C	X APPROXIMATION ALLOWS THIS SKIP
C
	IF(BS(J1) .NE. BS(K1)) GO TO 270
	XL=YL(J1)
	A(I,J)=(XL-Z(I1,J1))*(XL+Z(I1,K1))*(XL+Z(I1,J1)+1.)
	1 *(XL-Z(I1,K1)+1.)
	A(I,J)=0.5*YJ(J1,K1)*SQRT(A(I,J))
	A(J,I)=A(I,J)
270	CONTINUE
280	CONTINUE
C
C	CALL MATRIX DIAGONALIZATION (JACOBI METHOD, BUT ANY OTHER
C	 SUBROUTINE MAY BE SUBSTITUTED)
C
	CALL EIGEN2(N,R8)
	IF(H .EQ. L9) GO TO 360
C
C	COMPUTE INTENSITIES AND ENERGIES, WRITE THEM TO DISK FILE
C	(THIS SECTION SKIPPED FOR H=L9)
C
	DO 350 I=1,N6
	DO 340 J=1,N
	I7=0.
	DO 310 K=1,N
310	I7=I7+T(I,K)*S(K,J)
	IF(ABS(I7) .LT. I2) GO TO 340
	N7=N7+1
	I7=I7*I7
	I4=I4+I7
	E7=R(I)-A(J,J)
	IF(E7 .GE. E1) GO TO 320
	E1=E7
	GO TO 330
320	IF(E7 .LE. E2) GO TO 330
	E2=E7
330	WRITE(IO3,600) E7,I7
340	CONTINUE
350	CONTINUE
C
C	TEST IF ESCAPE KEY STRUCK DURING DIAGONALIZATION
C
360	N6=N
	I6=I0
	ESCAP=0
	IF(H .LE. -L9) RETURN
	H=H-1.
	CALL TTYTST (ESCAP,KDUM)
	IF(ESCAP .EQ. 32) RETURN
	GO TO 100
C
600	FORMAT(2E15.8)
	END

	SUBROUTINE ORDER (N7,N9,R0,KWIT)
C
C
C	READS ENERGIES AND INTENSITIES FROM DISK FILE, SORTS THEM,
C	 DEGRADES RESOLUTION IF NEEDED TO COMPRESS TO NO MORE THAN
C	 N9 LINES
C
	COMMON /SPEC/ E(500),I(500),X(500),Y(0/510)
C
	COMMON /IODEV/ IO2,IO3,IO4,IO5,IO6
C
	REAL I,I7
	INTEGER X
C
C	N7	COUNTER FOR NUMBER OF LINES
C	N9	STORAGE LIMIT
C	E7	ENERGY FOR LINE
C	I7	INTENSITY FOR LINE
C
C	FNR	ROUNDING FUNCTION TO COMPRESS SPECTRUM
C
	FNR(XX)=AINT(XX*R0+SIGN(0.5,XX))/R0
C
C	READ NEW LINE, INSERT INTO PROPER SEQUENCE
C
	KWIT=1
	N7=1
	READ(IO2,500,END=300) E7,I(1)
	E(1)=FNR(E7)
C
100	READ(IO2,500,END=300) E7,I7
	E7=FNR(E7)
	DO 120 L=1,N7
	IF(E7-E(L)) 110,140,120
110	T=E(L)
	E(L)=E7
	E7=T
	T=I(L)
	I(L)=I7
	I7=T
120	CONTINUE
C
130	IF(N7 .GE. N9) GO TO 200
	N7=N7+1
	E(N7)=E7
	I(N7)=I7
	GO TO 100
140	I(L)=I(L)+I7
	GO TO 100
C
C	DEGRADE RESOLUTION IF OUT OF STORAGE
C
200	KWIT=-1
	R0=0.5*R0
	E(1)=FNR(E(1))
	E7=FNR(E7)
	K=1
	DO 220 L=2,N7
	E(L)=FNR(E(L))
	IF(E(L) .NE. E(K)) GO TO 210
	I(K)=I(K)+I(L)
	GO TO 220
210	K=K+1
	E(K)=E(L)
	I(K)=I(L)
220	CONTINUE
	N7=K
	IF(E7 .EQ. E(N7)) GO TO 140
	GO TO 130
C
300	REWIND IO2
	RETURN
C
500	FORMAT(2E15.8)
	END

	SUBROUTINE SHAPE (N7,GP,H1,T1,T2,Y9)
C
C
C	USES LORENTZ LINESHAPE FUNCTION TO CONVERT LINE SPECTRUM
C	 FROM E(N7) TO E(1) INTO SIMULATED SPECTRUM
C
	COMMON /SPEC/ E(500),I(500),X(500),Y(0/510)
C
	COMMON /PLOT1/ KWIT,L,X0,X1,X2,X3,X4,X5,Y0,Y1,Y2,Y3,Y4,Y5
C
	INTEGER X,X2,X3,X4,X5,Y2,Y3,Y4,Y5
	REAL I,LO
C
C	X6	NUMBER OF BITS FOR 1 HZ
C	B4	CUT-OFF CRITERION FOR TAILS OF PEAKS
C	LO	LORENTZ AMPLITUDE FACTOR
C
C	INDEX LINES DIGITALLY ALONG X AXIS
C
	X6=FLOAT(X5)/(X1-X0+1.0E-5)
	DO 110 J=1,N7
	XT=X6*(E(J)-X0)
110	X(J)=IFIX(XT+SIGN(0.5,XT))
C
C	ZERO PREVIOUS PLOT
C
	DO 120 J=0,X5
120	Y(J)=0.
C
C	CONSTANTS FOR LORENTZ FUNCTION
C
	B1=GP*H1*T2
	B2=1.+GP*GP*H1*H1*T1*T2
	B3=(T2*T2)/(X6*X6)
	LO=B1/B2
	J1=0
	Y9=0.
C
C	ADJUST HEIGHTS AT PEAK CENTERS
C
	DO 130 J=1,N7
	MT=X(J)
	IF((MT .LT. 0) .OR. (MT .GT. X5)) GO TO 130
	YY=Y(MT)+LO*I(J)
	IF(YY .GT. Y9) Y9=YY
	Y(MT)=YY
130	CONTINUE
C
	B4=1.0E-5 + 2.*LO/FLOAT(Y5)
140	J1=J1+1
	LO=B1/(B2+B3*FLOAT(J1)*FLOAT(J1))
	IF(LO .LT. B4) RETURN
C
C	ADD OFF-CENTER AMPLITUDES FOR EVERY LINE
C	 (NOTE THAT J1 ALTERNATES LEFT AND RIGHT SIDES OF PEAK)
C
	DO 160 J=1,N7
	MT=X(J)
	DO 160 K=1,2
	K1=MT+J1
	IF((K1 .LT. 0) .OR. (K1 .GT. X5)) GO TO 150
	YY=Y(K1)+LO*I(J)
	IF(YY .GT. Y9) Y9=YY
	Y(K1)=YY
150	J1=-J1
160	CONTINUE
	GO TO 140
	END

	SUBROUTINE YPLOT (Y,P44)
C
	IMPLICIT INTEGER (A-Z)
C
C	DRAWS A LINE FROM PREVIOUS POINT TO THIS ONE AT SPEED P44
C
C	CALLING PROGRAM SETS L,X0,X1,X2,X5,Y0,Y1,Y5,
C	 BUT IT MUST NOT CHANGE OTHER VARIABLES IN BLOCK PLOT1
C
	COMMON /PLOT1/ KWIT,L,X0,X1,X2,X3,X4,X5,Y0,Y1,Y2,Y3,Y4,Y5
C
C	P	BUFFER CONTAINING CHARACTER CODES FOR PLOT COMMANDS
C	P0	NUMBER OF SMALL MOVEMENTS TO NEW POINT
C	P1	INDEX
C	P2	INDEX
C	P3	CHARACTER CODE (DECIMAL FOR 7-BIT ASCII)
C	P4(4)	NUMBER OF BITS PER SMALL MOVEMENT
C	P9	NUMBER OF BUFFER CELLS FILLED
C
	COMMON /PLOT2/ P(-2/66),P0,P1,P2,P3,P4,P9
C
	REAL X0,X1,Y,Y0,Y1
C
C	CONVERT Y TO INTEGER REPRESENTATION Y0 .LE. Y .LE. Y1
C
	P4=P44
	Y2=IFIX(FLOAT(Y5)*(Y-Y0)/(Y1-Y0+1.0E-5))
	IF(Y2 .LE. 0) Y2=0
	IF(Y2 .GT. Y5) Y2=Y5
C
C	COMPUTE DISPLACEMENT, ADJUST PEN VARIABLE (0=DOWN)
C
	X3=X2-X3
	Y3=Y2-Y3
	IF(L .NE. 0) L=64
C
C	VECTOR GENERATOR, P4=MAX MOVEMENT PER STEP
C
	P0=1+IABS(X3)/P4+IABS(Y3)/P4
	DO 100 P1=1,P0
	P9=P9+3
	X4=X2-(P0-P1)*X3/P0
	Y4=Y2-(P0-P1)*Y3/P0
	P3=64+X4/8
	IF(P3 .GE. 127) P3=63
	P(P9-2)=P3
	P3=64+Y4/8
	IF(P3 .GE. 127) P3=63
	P(P9-1)=P3
	P3=L+Y4-8*(Y4/8)+8*(X4-8*(X4/8))
	IF(P3 .GE. 127) P3=126
	P(P9)=P3
	IF(P9 .GE. 66) CALL PCTRL(2)
	IF(KWIT .EQ. 32) RETURN
100	CONTINUE
C
C	SET X3,Y3 = DESTINATION POINT, DROP PEN
C
	X3=X2
	Y3=Y2
	L=0
	RETURN
	END

	SUBROUTINE PCTRL(N)
C
	IMPLICIT INTEGER (A-Z)
C
C	N=1, ENTER PLOT MODE, TURN OFF TTY ECHO
C	N=2, PLOT CONTENTS OF BUFFER, TEST IF ESCAPE KEY STRUCK
C	N=3, RETURN TO PRINT MODE, RESTORE TTY ECHO
C
	COMMON /PLOT1/ KWIT,L,X0,X1,X2,X3,X4,X5,Y0,Y1,Y2,Y3,Y4,Y5
C
	COMMON /PLOT2/ P(-2/66),P0,P1,P2,P3,P4,P9
C
	REAL X0,X1,Y0,Y1
C
	GO TO (10,20,20) N
C
10	P(-2)=13
	P(-1)=10
	P(0)=16
	P9=0
	CALL NECHO
	RETURN
C
20	KWIT=0
	CALL TTYTST(KWIT,K)
	IF(KWIT .EQ. 32) GO TO 24
	DO 22 P2=-2,P9
22	CALL CHOUT(P(P2))
24	P9=0
	IF((N .EQ. 2) .AND. (KWIT .NE. 32)) RETURN
C
	P(0)=0
	DO 28 P2=-2,0
28	CALL CHOUT(P(P2))
	CALL ECHO
	RETURN
	END

	SUBROUTINE EIGEN2(N,R8)
C
C
C	MATRIX DIAGONALIZATION BY METHOD OF JACOBI(1846) USING
C	PIVOT SEARCH TECHNIQUE OF F. J. CORBATO(J.ASSOC. COMPUT.
C	MACH. 10, 123-5 (1963)) BUT NOTATION OF J. GREENSTADT
C	(PP. 84-91 IN "MATHEMATICAL METHODS FOR DIGITAL COMPUTERS"
C	ED. BY A. RALSTON & H. S. WILF (WILEY, 1960))
C
	FNA(A1,A2)=C1*(A1-A2*T1)
	FNB(A1,A2)=C1*(A2+A1*T1)
C
	COMMON /CALC/ A(35,35),S(35,35),T(35,35),R(35),X(35),MY(35),
	1 OP(128),Z(128,0/7),GZ(7),YJ(7,7),YL(7),BS(0/7),W(15),WV(15)
C
	INTEGER P,Q
C
C	SET S TO IDENTITY MATRIX
C	FIND BIGGEST ELEMENT IN EACH COLUMN OF A
C
	S(1,1)=1.0
	IF(N.LT.2) GO TO 2760
	DO 2640 J=2,N
	CALL BIGJ(J)
	S(J,J)=1.0
	DO 2630 I=1,J-1
	S(I,J)=0.0
	S(J,I)=0.0
2630	CONTINUE
2640	CONTINUE
C
C	FIND PIVOT
C
2660	V2=X(2)
	P=MY(2)
	Q=2
	IF(N .LT. 3) GO TO 2750
	DO 2740 J=3,N
	IF(V2.GE.X(J)) GO TO 2740
	V2=X(J)
	P=MY(J)
	Q=J
2740	CONTINUE
C
C	TEST FOR CONVERGENCE
C
2750	IF(V2.GT.R8) GO TO 2770
2760	RETURN
C
C	CALCULATE TRIG FUNCTIONS
C
2770	V1=A(P,P)
	V2=A(P,Q)
	V3=A(Q,Q)
	V4=V1-V3
	T1=-2*V2*SIGN(1.,V4)/(ABS(V4)+SQRT(V4*V4+4*V2*V2))
	IF(V4.EQ.0.) T1=-1.0
	T2=T1*T1
	C2=1/(1+T2)
	C1=SQRT(C2)
C
C	ROTATE MATRIX; REVISE X & MY AS NECESSARY
C
	V5=2*T1*V2
	A(P,P)=C2*(V1-V5+T2*V3)
	A(Q,Q)=C2*(V3+V5+T2*V1)
	A(P,Q)=0.0
C
	IF(Q.LT.(P+2)) GO TO 3010
	DO 3000 J=P+1,Q-1
	IF(MY(J).NE.P) GO TO 3000
	A1=A(P,J)
	A(P,J)=0.0
	CALL BIGJ(J)
	A(P,J)=A1
3000	CONTINUE
C
3010	IF((Q+1).GT.N) GO TO 3100
	DO 3090 J=Q+1,N
	IF((MY(J) .NE. P) .AND. (MY(J) .NE. Q)) GO TO 3090
	A1=A(P,J)
	A2=A(Q,J)
	A(P,J)=0.0
	A(Q,J)=0.0
	CALL BIGJ(J)
	A(P,J)=A1
	A(Q,J)=A2
3090	CONTINUE
C
3100	IF(P.LT.2) GO TO 3160
	DO 3150 I=1,P-1
	A1=A(I,P)
	A(I,P)=FNA(A1,A(I,Q))
	A(I,Q)=FNB(A1,A(I,Q))
3150	CONTINUE
C
3160	IF(Q.LT.(P+2)) GO TO 3240
	DO 3230 J=P+1,Q-1
	A1=A(P,J)
	A(P,J)=FNA(A1,A(J,Q))
	A(J,Q)=FNB(A1,A(J,Q))
	IF(X(J).GE.ABS(A(P,J))) GO TO 3230
	MY(J)=P
	X(J)=ABS(A(P,J))
3230	CONTINUE
C
3240	IF((Q+1).GT.N) GO TO 3340
	DO 3330 J=Q+1,N
	A1=A(P,J)
	A(P,J)=FNA(A1,A(Q,J))
	A(Q,J)=FNB(A1,A(Q,J))
	IF(X(J).GE.ABS(A(P,J))) GO TO 3320
	MY(J)=P
	X(J)=ABS(A(P,J))
3320	IF(X(J).GE.ABS(A(Q,J))) GO TO 3330
	MY(J)=Q
	X(J)=ABS(A(Q,J))
3330	CONTINUE
C
3340	CALL BIGJ(P)
	CALL BIGJ(Q)
C
C	ROTATE EIGENVECTOR MATRIX
C
	DO 3440 I=1,N
	A1=S(I,P)
	S(I,P)=FNA(A1,S(I,Q))
	S(I,Q)=FNB(A1,S(I,Q))
3440	CONTINUE
	GO TO 2660
	END

	SUBROUTINE BIGJ(J)
C
C
C	FINDS BIGGEST ELEMENT IN COLUMN J (ABOVE DIAGONAL)
C	SAVES MAGNITUDE IN X(J), ROW INDEX IN MY(J)
C
	COMMON /CALC/ A(35,35),S(35,35),T(35,35),R(35),X(35),MY(35),
	1 P(128),Z(128,0/7),GZ(7),YJ(7,7),YL(7),BS(0/7),W(15),WV(15)
C
	X(J)=ABS(A(1,J))
	MY(J)=1
	IF(J .LT. 3) RETURN
	DO 3510 I=2,J-1
	IF(X(J).GE.ABS(A(I,J))) GO TO 3510
	MY(J)=I
	X(J)=ABS(A(I,J))
3510	CONTINUE
	RETURN
	END