Program RegressGT * A regression program to fit Gastherm logk's v. T to an equation of the form: * LOG K (T(K) = A + B/T + CT + D/T**2 + E(LOG(T)) * This program reads the log k data from a Gastherm data file and creates * a new Gastherm data file containing the regression coefficients * Written by M.H. Reed, 1976, modified by N.Spycher 1983. * Modified by Bob Symonds, April 1986 to read Gastherm, and April 1990 * to use the equation above, replacing the polynomial form. * Modified for readability and fewer compile warnings, J.Palandri 2023. * Only the main program has been modified significantly. Note however * that more meaningful variable names have been propagated through * the subroutines. implicit none real*8 T_K, Log_K, TK(20), LgK(20), TempK(20), LogK(20) integer iloop, jloop, inp1, iout1, iout2 integer n, nd, ne, jk, l1, l2 character*8 stoch(12), name character*36 ref inp1 = 6 !input data iout1 = 7 !output coefficient file iout2 = 8 !output regression analysis open (inp1, file = 'RegrunGT.dat', status = 'old') open (iout1, file = 'RegoutGT.dat', status = 'unknown') open (iout2, file = 'RegpltGT.dat', status = 'unknown') nd = 13 !# of logK's in each record n = 4 !rxn order always 4 here (was variable for Soltherm polynomial) read (inp1,'(i5)') ne !# of input file records at the top of inp1 * Create values TempK(i), the temperatures at which log(k)'s are tabulated TempK(1) = 298.15 TempK(2) = 373.15 jk = 2 do 3 iloop = 3,nd TempK(iloop) = TempK(jk) + 100. jk = iloop 3 continue do 6 iloop = 1,ne read (inp1,'(12a8)') (stoch(jloop), jloop=1,12) read (inp1,'(a8,8f8.3,a36)') name, (LogK(jloop),jloop=1,8), ref read (inp1,'(8x,8f8.3)') (LogK(jloop),jloop=9,nd) read (inp1,*) * l1,l2 are for options for calculating regression error parameters l1 = 0 l2 = 1 do 4 jloop = 1,nd T_K = TempK(jloop) Log_K = LogK(jloop) TK(jloop) = T_K LgK(jloop) = Log_K 4 continue call least (n,nd,l1,TK,LgK,l2,name,stoch,LogK,ref) 6 continue end *===================================================================== SUBROUTINE LEAST(N,ND,L1,TK,LGK,L2,NAME,STOCH,LOGK,REF) *===================================================================== CHARACTER*8 NAME,STOCH(12) character*36 REF REAL*8 EX(5,13), TK(20), LGK(20), S1, S2, S3, S4, S5, S6 REAL*8 XBAR, YBAR, FM, YDIF, FN REAL*8 A(20,21), C(20), D(20), LOGK(20) IOUT1 = 7 IOUT2 = 8 IF(L2)201,201,200 200 WRITE (IOUT2,5) NAME 5 FORMAT(////,2X,'LOGK''S FOR: ',A8 , + ' REGRESSION TO FIT LOGK(T(K)) = ', + 'A + B/T + CT + D/T**2 + E(LOG(T)' ,/104('*')) 201 FN=N FM=ND NP1=N+1 * 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+TK(I) S2=S2+LGK(I) 100 CONTINUE XBAR=S1/FM YBAR=S2/FM IF(L1)101,103,103 101 DO 102 I=1,ND TK(I)=TK(I)-XBAR LGK(I)=LGK(I)-YBAR 102 CONTINUE 103 CONTINUE * EX1,EX2,EX3,EX4 ARE THE EQUATIONS USED IN THE REGRESSION DO 104 J=1,ND EX(1,J)=1. EX(2,J)=1./TK(J) EX(3,J)=TK(J) EX(4,J)=1./(TK(J)*TK(J)) 104 EX(5,J)=DLOG10(TK(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 * RESET A(1,1) *109 A(1,1)=FM S1=0.0 DO 110 I=1,ND 110 S1=S1+LgK(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)*LgK(J) 112 A(I,N+2)=S1 I=1 J=NP1 CALL SLNQ(A,J,I,D) IF(L2)203,203,202 202 WRITE (IOUT2,10) (D(I),I=1,NP1) WRITE (IOUT1,14) (STOCH(I),I=1,12) WRITE (IOUT1,15) NAME,(LOGK(I),I=1,8),REF WRITE (IOUT1,15) NAME, (LOGK(I),I=9,ND) WRITE (IOUT1,16) NAME, (D(I), I=1,NP1) WRITE (IOUT2,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=LgK(I)-S1 write (iout2,12) TK(I),LgK(I),S1,YDIF S2=S2+YDIF*YDIF IF(L1)115,116,116 115 S3=S3+TK(I)*LgK(I) S4=S4+TK(I)*TK(I) S5=S5+LgK(I)*LgK(I) GO TO 117 116 S3=S3+(TK(I)-XBAR)*(LgK(I)-YBAR) S4=S4+(TK(I)-XBAR)*(TK(I)-XBAR) S5=S5+(LgK(I)-YBAR)*(LgK(I)-YBAR) 117 DO 119 J=1,NP1 IF(TK(I))118,119,118 118 C(J)=C(J)+(YDIF/(TK(I)**J))*(YDIF/(TK(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 (iout2,13) S2,S3,S6,(C(I),I=1,NP1) * PLOT INITIALISATION * M=1 L=31 CALL PLOT(D,L,M) 10 FORMAT(/1X,12HCOEFFICIENTS//(1X,8E16.8)) 11 FORMAT (/12X,'TEMP IN DEG K',17X,'EXP LOG K',15X,'CALC LOG K', +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 +NT,E16.8,5X,11HCHI SQUARED,E16.8//1X,31HSTANDARD ERRORS IN COEFFIC +IENTS//(1X,8E16.8)) 14 FORMAT(12A8) 15 FORMAT(A8,8F8.3,A36) 16 FORMAT (A8,5E14.7) 203 RETURN END *===================================================================== SUBROUTINE SLNQ(A,N,J,TK) *===================================================================== * SOLUTION OF N LINEAR EQ. WITH N AN. AND DET. EV. DOUBLE PRECISION A(20,21),TK(20),Y,B,C,D IOUT2 = 8 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 TK(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)*TK(L) 112 CONTINUE TK(K)=A(K,N1)-B K=K-1 IF(K)114,113,111 113 RETURN 114 IF(J)115,116,115 115 write (iout2,10) STOP 116 Y=0.0 RETURN END *===================================================================== SUBROUTINE PLOT(D,L,M) *===================================================================== * AUTHOR K.J. JOHNSON INTEGER CHAR(4) DOUBLE PRECISION TK(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./ IOUT2 = 8 * CALCULATE Y FOR X = 25 TO 1227 STEP 40 TK(1) = 300. DO 300 I=1,L XX=TK(I) Y(I,M)=D(1)+D(2)/XX+D(3)*XX+D(4)/(XX*XX)+D(5)*DLOG10(XX) TK(I+1)=TK(I)+40. 300 CONTINUE * 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 (iout2,11) YMIN,YMID,YMAX * PLOT ONE LINE AT A TIME 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 (iout2,22) TK(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