C REGMK IS A SPECIAL POLYNOMIAL REGRESSION PROGRAM TO FIT LOG K'S C CALCULATED FROM MAIER-KELLY HEAT CAPICITY EQUATIONS TO AN C EQUATION OF THE FORM: C LOG K (T(K) = A + B/T + CT + D/T**2 + E(LOG(T)). C THE PROGRAM READS THE LOG K DATA FROM THE GASTHERM DATA FILE AND C CREATES A NEW GASTHERM DATA FILE CONTAINING THE REGRESSION C COEFFICIENTS. C WRITTEN BY M.H. REED, 1976, MODIFIED BY N.SPYCHER 1983 C MODIFIED BY BOB SYMONDS, APRIL 1986 TO READ GASTHERM C BOB SYMONDS (4,90) TO USE ABOVE EQUATION C DEVICE NUMBERS: 10 INPUT DATA C 6 OUTPUT REGRESSION ANALYSIS C 7 OUTPUT COEFFICIENT FILE C DOUBLE PRECISION NAME,STOCH,TEXT,ETOIL,REF,XXX,YYY DIMENSION STOCH(10) DOUBLE PRECISION X(20),Y(20),XQ(20), YQ(20) C NA IS THE TOTAL NUMBER OF Y VALUES, N IS REGRESSION ORDER C (ALWAYS 4TH ORDER FOR GASTHERM), ND IS # OF Y VALUES REGRESSED NA = 13 N = 4 DATA ETOIL/'HF(G) '/ C C OPEN STATEMENTS FOR IBM DOS MACHINES C OPEN (UNIT = 10, FILE = 'REGIN.DAT', STATUS = 'OLD') OPEN (UNIT = 6, FILE = 'REGMK.DAT', STATUS = 'UNKNOWN') OPEN (UNIT = 7, FILE = 'REGOUT.DAT', STATUS = 'UNKNOWN') C READ (10,100) NE 2 READ (7,105) TEXT IF(TEXT.NE.ETOIL) GO TO 2 C READ X VALUES (THE TEMPERATURES FOR WHICH LOG(K)'S ARE GIVEN XQ(1) = 298.15 XQ(2) = 373.15 JK = 2 DO 3 I = 3,NA XQ(I) = XQ(JK) + 100. JK = I 3 CONTINUE C DO 6 NZ = 1,NE READ(10,110) (STOCH(I), I=1,10) READ (10,115) NAME, (YQ(K),K=1,8), REF READ (10,116) (YQ(K),K=9,NA) READ (10,116) C L1,L2 ARE FOR OPTIONS FOR CALCULATING REGRESSION ERROR PARAMETERS L1= 0 L2=1 ND = 13 DO 4 KZ = 1,ND XXX = XQ(KZ) YYY = YQ(KZ) X(KZ) = XXX Y(KZ) = YYY 4 CONTINUE CALL LEAST (N,ND,L1,X,Y,L2,NAME,STOCH,YQ,NA,REF) 6 CONTINUE 100 FORMAT(I5) 105 FORMAT(A8,72X) 110 FORMAT(10A8) 115 FORMAT (A8,8F8.3,A8) 116 FORMAT (8X,8F8.3) 500 STOP END SUBROUTINE LEAST(N,ND,L1,X,Y,L2,NAME,STOCH,YQ,NA,REF) DOUBLE PRECISION NAME,STOCH,REF DOUBLE PRECISION EX(5,13),X(20),Y(20),S1,S2,S3,S4,S5,S6 DOUBLE PRECISION XBAR,YBAR,FM,YDIF,FN DIMENSION STOCH(10) DOUBLE PRECISION A(20,21),C(20),D(20),YQ(20) 5 FORMAT(1H1,6('*'),'LOG K FOR : ',A8 ,' POLYNOMIAL REGRESSION 1WITH LOG K(T(K)) = A + B/T + CT + D/T**2 + E(LOG(T)',/1X,104('*')) 10 FORMAT(/1X,12HCOEFFICIENTS//(1X,8E16.8)) 11 FORMAT (/12X,'TEMP IN DEG K',17X,'EXP LOG K',15X,'CALC LOG K', 1 11X,'EXP-CALC LOG K'/) 12 FORMAT(1X,3(15X,F10.4),15X,E12.5) 13 FORMAT(/1X,18HSTANDARD DEVIATION,E16.8,5X,23HCORRELATION COEFFICIE 1NT,E16.8,5X,11HCHI SQUARED,E16.8//1X,31HSTANDARD ERRORS IN COEFFIC 2IENTS//(1X,8E16.8)) 14 FORMAT(10A8) 15 FORMAT(A8,8F8.3,A8) 16 FORMAT (A8,5E14.7) IF(L2)201,201,200 200 WRITE (6,5) NAME 201 FN=N FM=ND NP1=N+1 C INITIALIZE VARIABLES S1=0.0 S2=0.0 DO 90 I=1,20 DO 80 J=1,21 A(I,J) = 0.0 80 CONTINUE 90 CONTINUE DO 100 I=1,ND S1=S1+X(I) 100 S2=S2+Y(I) XBAR=S1/FM YBAR=S2/FM IF(L1)101,103,103 101 DO 102 I=1,ND X(I)=X(I)-XBAR 102 Y(I)=Y(I)-YBAR 103 CONTINUE C C EX1,EX2,EX3,EX4 ARE THE EQUATIONS USED IN THE REGRESSION C DO 104 J=1,ND EX(1,J)=1. EX(2,J)=1./X(J) EX(3,J)=X(J) EX(4,J)=1./(X(J)*X(J)) 104 EX(5,J)=DLOG10(X(J)) DO 106 I=1,NP1 DO 105 J=1,NP1 S1=0. DO 109 K=1,ND S1 = S1 + EX(I,K)*EX(J,K) 109 CONTINUE A(I,J) = S1 105 CONTINUE 106 CONTINUE C C RESET A(1,1) C C109 A(1,1)=FM C S1=0.0 DO 110 I=1,ND 110 S1=S1+Y(I) A(1,N+2)=S1 DO 112 I=2,NP1 S1=0.0 DO 111 J=1,ND 111 S1=S1+EX(I,J)*Y(J) 112 A(I,N+2)=S1 I=1 J=NP1 CALL SLNQ(A,J,I,D) IF(L2)203,203,202 202 WRITE (6,10) (D(I),I=1,NP1) WRITE(7,14) (STOCH(I),I=1,10) WRITE(7,15) NAME,(YQ(I),I=1,8),REF WRITE(7,15) NAME, (YQ(I),I=9,NA) WRITE (7,16) NAME, (D(I), I=1,NP1) WRITE (6,11) DO 113 I=1,NP1 113 C(I)=0.0 S2=0.0 S3=0.0 S4=0.0 S5=0.0 S6=0.0 DO 120 I=1,ND S1=0.0 S1=S1+D(1) DO 114 J=2,NP1 114 S1=S1+EX(J,I)*D(J) YDIF=Y(I)-S1 WRITE (6,12) X(I),Y(I),S1,YDIF S2=S2+YDIF*YDIF IF(L1)115,116,116 115 S3=S3+X(I)*Y(I) S4=S4+X(I)*X(I) S5=S5+Y(I)*Y(I) GO TO 117 116 S3=S3+(X(I)-XBAR)*(Y(I)-YBAR) S4=S4+(X(I)-XBAR)*(X(I)-XBAR) S5=S5+(Y(I)-YBAR)*(Y(I)-YBAR) 117 DO 119 J=1,NP1 IF(X(I))118,119,118 118 C(J)=C(J)+(YDIF/(X(I)**J))*(YDIF/(X(I)**J)) 119 CONTINUE S6=S6+(YDIF*YDIF)/S1 120 CONTINUE S2=DSQRT(S2/(FM-1.0)) S3=S3/(DSQRT(S4*S5)) DO 121 I=1,NP1 121 C(I)=DSQRT(C(I)/FN) WRITE (6,13) S2,S3,S6,(C(I),I=1,NP1) C PLOT INITIALISATION C M=1 L=31 CALL PLOT(D,L,M) 203 RETURN END SUBROUTINE SLNQ(A,N,J,X) C SOLUTION OF N LINEAR EQ. WITH N AN. AND DET. EV. DOUBLE PRECISION A(20,21),X(20),Y,B,C,D 10 FORMAT(1X,27HSOLUTION OF LNEQ NOT UNIQUE) 100 Y=1.0 N1=N+1 N2=N-1 DO 109 I=1,N2 M=I B=A(I,I) DO 102 L=I,N C= DABS(B) IF(C- DABS (A(L,I)))101,102,102 101 M=L B=A(L,I) 102 CONTINUE IF(B)103,114,103 103 CONTINUE IF(I-M)104,105,104 104 Y=-Y 105 DO 106 L=I,N1 B=A(I,L) A(I,L)=A(M,L) A(M,L)=B 106 CONTINUE B=A(I,I) DO 107 L=I,N1 A(I,L)=A(I,L)/B 107 CONTINUE Y=Y*B K=I+1 DO 108 M=K,N D=A(M,I) DO 108 L=I,N1 A(M,L)=A(M,L)-D*A(I,L) 108 CONTINUE 109 CONTINUE Y=Y*A(N,N) IF(J)110,113,110 110 X(N)=A(N,N1)/A(N,N) K=N-1 111 M=K+1 B=0.0 DO 112 L=M,N B=B+A(K,L)*X(L) 112 CONTINUE X(K)=A(K,N1)-B K=K-1 IF(K)114,113,111 113 RETURN 114 IF(J)115,116,115 115 WRITE (6,10) STOP 116 Y=0.0 RETURN END SUBROUTINE PLOT(D,L,M) C C AUTHOR K.J. JOHNSON C INTEGER CHAR(4) DOUBLE PRECISION X(50),Y(50,3),D(20),XX,YMAX,YMID,YMIN,YSCALE DOUBLE PRECISION WIDTH DIMENSION LINE(51),NN(3) DATA LINE/51*' '/,CHAR/'+','-','*',' '/ DATA LENGTH,WIDTH/50,50./ C CALCULATE Y FOR X = 25 TO 1227 STEP 40 X(1) = 300. DO 300 I=1,L XX=X(I) Y(I,M)=D(1)+D(2)/XX+D(3)*XX+D(4)/(XX*XX)+D(5)*DLOG10(XX) X(I+1)=X(I)+40. 300 CONTINUE C FIND MINIMUM AND MAXIMUM Y AND SET SCALING FACTORS YMIN=Y(1,1) YMAX=Y(1,1) DO 10 I=1,L DO 10 J=1,M IF(Y(I,J).GT.YMAX)YMAX=Y(I,J) 10 IF (Y(I,J).LT.YMIN) YMIN=Y(I,J) YSCALE=WIDTH/(YMAX-YMIN) YMID=(YMAX+YMIN)/2. ISKIP=L/LENGTH+1 IF (MOD(L,LENGTH).EQ.0) ISKIP=ISKIP-1 WRITE (6,11) YMIN,YMID,YMAX C C PLOT ONE LINE AT A TIME C DO 40 I=1,L,ISKIP DO 30 J=1,M NN(J)=(Y(I,J)-YMIN)*YSCALE + 1.5 30 LINE(NN(J))=CHAR(J) NMAX=AMAX0(NN(1),NN(2),NN(3)) WRITE (6,22) X(I), (LINE(JJ),JJ=1,NMAX) DO 25 J=1,M 25 LINE(NN(J))=CHAR(4) 40 CONTINUE 11 FORMAT (/3X,'TEMP',25X,' LOG (K) ',//,9X,G12.3,12X,G12.3,12X 1G12.3,/,14X,':',24X,':',24X,':',/, 213X,':+----+----+----+----+----+----+----+----+----+----+') 22 FORMAT (6X,F7.2,' :',51A1) RETURN END