PROGRAM POLY C C ***************************************************** C * POLY.ODL as follows * C * * C * * C .ROOT POLY/LB:POLY-FLIB-*(S1,S2,S3) * C S1: .FCTR POLY/LB:INDATA-FLIB * C S2: .FCTR POLY/LB:FIT:OUT-FLIB * C S3: .FCTR POLY/LB:GRAPH-PLIB-MLIB-FLIB * C FLIB: .FCTR LB:[1,1]FOROTS/LB * C PLIB: .FCTR LB:[1,1]PLOT/LB * C MLIB: .FCTR LB:[1,1]MURDLIB/LB * C .END * C ***************************************************** C C C LU5 - CON: FOR INTERACTION OUTPUT (NULL: FOR BATCH MODE) C LU6 - CON: FOR INTERACTION INPUT (CARD: FOR BATCH MODE) C LU1 - INPUT OF DATA C LU2 - OUTPUT OF ALPHA C LU8 - INTERMEDIATE PLOT FILE (IPF) C LOGICAL SCRN INTEGER ANS BYTE INF(20),OUTF(20),IDEV C COMMON IBUF(500),Z(502) COMMON /A/C1,E1,F1,D1,YE(502),IORDA,I2,Y(202),A(11),X(202) COMMON /B/ X1,X2,X3,X4,X5,X6,Y1,Y2,Y3,Y4 COMMON /C/ ALL1,ALL2,ALL7,ALL8,ALLXY,ALXY,SCRN C EQUIVALENCE (A(1),A1),(A(2),A2),(A(3),A3),(A(4),A4),(A(5),A5), 1 (A(6),A6),(A(7),A7),(A(8),A8),(A(9),A9),(A(10),A10), 1 (A(11),A11) C DATA LU5/5/,LU6/6/ DATA SCRN/.FALSE./ DATA INF/'T','I',':',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ', 1 ' ',' ',' ',' ',' ',' ',' '/ C C OPEN FILES C PRINT 4300 4300 FORMAT('$Data input from screen (Y/N) ?') ACCEPT 4400,IDEV 4400 FORMAT(A1) IF(IDEV.EQ.'N'.OR.IDEV.EQ.'n')GOTO 4500 SCRN=.TRUE. GOTO 4600 4500 PRINT 4000 4000 FORMAT('$Enter name of input file ->') ACCEPT 4100,INF 4100 FORMAT(20A1) 4600 INF(20)=0 OPEN(UNIT=1,NAME=INF,TYPE='OLD') PRINT 4200 4200 FORMAT('$Enter name of output device/file ->') ACCEPT 4100,OUTF OPEN(UNIT=2,NAME=OUTF) C C WRITE PROGRAM ANNOUNCEMENT C WRITE (LU5,301) 301 FORMAT (43H0 THE CURVE-FITTING PROGRAM IS NOW RUNNING. ) CALL INDATA 322 CALL FIT C 321 WRITE (LU5,320) 320 FORMAT (' DO YOU WANT A GRAPHICS PLOT? 1=YES, 2=NO.') READ (LU6,*) ANS IF (ANS.EQ.2) GO TO 323 IF (ANS.NE.1) GO TO 321 CALL GRAPH 323 WRITE (LU5,324) 324 FORMAT (' DO YOU WANT ANOTHER FIT ON THE DATA? 1=YES, 2=NO.') READ (LU6,*) ANS IF (ANS.EQ.2) GO TO 326 IF (ANS.NE.1) GO TO 323 C C ANOTHER FIT ON THE SAME DATA. C 399 WRITE (LU5,400) 400 FORMAT (7H0ENTER:/ 1 28H 1 FOR LINEAR THROUGH ORIGIN / 2 28H 2 FOR LINEAR WITH INTERCEPT / 3 31H 3 FOR QUADRATIC THROUGH ORIGIN / 4 31H 4 FOR QUADRATIC WITH INTERCEPT / 5 27H 5 FOR CUBIC THROUGH ORIGIN / 6 27H 6 FOR CUBIC WITH INTERCEPT / 7 31H 7 FOR EXPONENTIAL (LOG-LINEAR) / 8 22H 8 FOR POWER (LOG-LOG) ) C READ (LU6,*) IORDA IF (IORDA.LT.1.OR.IORDA.GT.8) GO TO 399 C GO TO 322 326 WRITE (LU5,325) 325 FORMAT (12H END OF RUN. ) STOP END SUBROUTINE INDATA C LOGICAL SCRN COMMON INBUF(500),Z(502) COMMON /A/ C1,E1,F1,D1,YE(502),IORDA,I2,Y(202),A(11),X(202) COMMON /B/ X1,X2,X3,X4,X5,X6,Y1,Y2,Y3,Y4 COMMON /C/ ALL1,ALL2,ALL7,ALL8,ALLXY,ALXY,SCRN C EQUIVALENCE (A(1),A1),(A(2),A2),(A(3),A3),(A(4),A4),(A(5),A5), 1 (A(6),A6),(A(7),A7),(A(8),A8),(A(9),A9),(A(10),A10), 2 (A(11),A11) C DATA AOUT /-99./ DATA LU5/5/,LU6/6/ C DO 2 I=1,11 2 A(I)=0.0 ALL1 =0. ALL2 =0. ALL7 =0. ALL8 =0. ALLXY=0. ALXY =0. C 21 WRITE(LU5,203) 203 FORMAT (7H0ENTER: / 1 28H 1 FOR LINEAR THROUGH ORIGIN / 2 28H 2 FOR LINEAR WITH INTERCEPT / 3 31H 3 FOR QUADRATIC THROUGH ORIGIN / 4 31H 4 FOR QUADRATIC WITH INTERCEPT / 5 27H 5 FOR CUBIC THROUGH ORIGIN / 6 27H 6 FOR CUBIC WITH INTERCEPT / 7 31H 7 FOR EXPONENTIAL (LOG-LINEAR) / 8 22H 8 FOR POWER (LOG-LOG) ) READ (LU6,*) IORDA IF (IORDA.LT.1.OR.IORDA.GT.8) GO TO 21 IF(.NOT.SCRN)GOTO 4000 WRITE (LU5,302) AOUT,AOUT 302 FORMAT (52H ENTER THE NUMBERS ONE X,Y PAIR PER LINE FREEFIELD. / 1 13H ENTER VALUES,F5.0,1H,,F5.0,24H TO SIGNIFY END OF DATA. ) C 4000 DO 3 I=1,202 READ (1,*,END=4) X(I),Y(I) IF (X(I).EQ.AOUT.AND.Y(I).EQ.AOUT) GO TO 4 3 CONTINUE C TOO MANY DATA POINTS WRITE (LU5,32) 32 FORMAT (49H0TOO MANY DATA POINTS. ONLY 1ST 200 WILL BE USED. ) I=201 4 I2=I-1 C**** ASSIGN MATRIX FACTORS DO 5 I=1,I2 X1=X(I) X2=X1*X1 X3=X2*X1 X4=X2*X2 X5=X2*X3 X6=X2*X4 Y1=Y(I) Y2=Y1*X1 Y3=Y1*X2 Y4=Y1*X3 A1=A1+X1 A2=A2+X2 A3=A3+X3 A4=A4+X4 A5=A5+X5 A6=A6+X6 A7=A7+Y1 A8=A8+Y2 A9=A9+Y3 A10=A10+Y4 A11=A11+1 C LOG VALUES ALLX =ALOG(X(I)) ALLY =ALOG(Y(I)) ALL1 =ALL1+ALLX ALL2 =ALL2+ALLX*ALLX ALL7 =ALL7+ALLY ALL8 =ALL8+ALLX*ALLY ALLXY=ALLXY+ALLX*ALLY ALXY =ALXY+X1*ALLY 5 CONTINUE C RETURN END SUBROUTINE FIT C LOGICAL SCRN COMMON IBUF(500),Z(502) COMMON /A/ C1,E1,F1,D1,YE(502),IORDA,I2,Y(202),A(11),X(202) COMMON /B/ X1,X2,X3,X4,X5,X6,Y1,Y2,Y3,Y4 COMMON /C/ ALL1,ALL2,ALL7,ALL8,ALLXY,ALXY,SCRN C EQUIVALENCE (A(1),A1),(A(2),A2),(A(3),A3),(A(4),A4),(A(5),A5), 1 (A(6),A6),(A(7),A7),(A(8),A8),(A(9),A9),(A(10),A10), 2 (A(11),A11) C DATA LU5/5/,LU6/6/ C C1 =0. D1 =0. E1 =0. F1 =0. C GO TO (10,8,9,7,12,13,14,16),IORDA C**** STRAIGHT LINE THRU ZERO 10 E1=A8/A2 GO TO 15 C**** STRAIGHT LINE ON INTERCEPT 8 E1=(A11*A8-A1*A7)/(A11*A2-A1*A1) F1=(A7*A2-A1*A8)/(A11*A2-A1*A1) GO TO 15 C**** QUADRATIC THRU ZERO 9 E1=(A9*A3-A8*A4)/(A3*A3-A2*A4) D1=(A8-E1*A2)/A3 GO TO 15 C**** QUADRATIC ON INTERCEPT 7 Z1=A7*A1-A11*A8 Z2=A1*A2-A11*A3 Z3=A2*A7-A11*A9 Z4=A1*A1-A11*A2 Z5=A2*A2-A11*A4 D1=(Z1*Z2-Z3*Z4)/(Z2*Z2-Z4*Z5) E1=(Z1-D1*Z2)/Z4 F1=(A7-E1*A1-D1*A2)/A11 GO TO 15 C**** CUBIC THRU ZERO 12 Z1=A4*A6-A5*A5 Z2=A3*A5-A4*A4 Z3=A4*A5-A3*A6 Z4=A2*A6-A4*A4 Z5=A3*A4-A2*A5 Z6=A2*A4-A3*A3 DEL=A2*Z1+A3*Z3+A4*Z2 IF (ABS(DEL).GT.1.E-5) GO TO 318 WRITE (LU5,316) RETURN 318 E1=(A8*Z1+A10*Z2+A9*Z3)/DEL D1=(A9*Z4+A8*Z3+A10*Z5)/DEL C1=(A10*Z6+A9*Z5+A8*Z2)/DEL GO TO 15 C**** CUBIC WITH INTERCEPT 13 Z1=A1*A4-A2*A3 Z4=A1*A5-A2*A4 Z6=A1*A3-A2*A2 Z7=A1*A6-A3*A4 Z8=A1*A5-A3*A3 Z9=A2*A5-A3*A4 Z10=A2*A3-A11*A5 Z11=A2*A2-A11*A4 Z12=A2*A6-A3*A5 Z13=A1*A2-A11*A3 Z14=A2*A4-A3*A3 Z15=A1*A1*A1 Z16=A1*A3-A11*A4 Z17=A1*A1-A11*A2 DEL=-A1*((A1**2-A11*A2)/Z15*(Z1*Z7-Z8*Z4) 1+(A2*A1-A11*A3)/Z15*(Z1*Z4-Z6*Z7) 1+(A3*A1-A11*A4)/Z15*(Z6*Z8-Z1*Z1)) IF (ABS(DEL).GT.1.E-5) GO TO 315 WRITE (LU5,316) 316 FORMAT (54H THE DETERMINANT IS TOO SMALL. THIS FIT WON'T BE DONE.) RETURN 315 F1=-A1*((A8*A1-A2*A7)/Z15*(Z1*Z7-Z8*Z4) 1+(A9*A1-A7*A3)/Z15*(Z1*Z4-Z6*Z7) 1+(A10*A1-A7*A4)/Z15*(Z6*Z8-Z1*Z1))/DEL E1=A2*((A8*A2-A7*A3)/A2**3*(Z9*Z10-Z11*Z12) 1+(A9*A2-A7*A4)/A2**3*(Z13*Z12-Z14*Z10) 1+(A10*A2-A7*A5)/A2**3*(Z14*Z11-Z13*Z9))/DEL D1=-A1*((A8*A1-A7*A2)/Z15*(Z8*Z16-Z13*Z7) 1+(A9*A1-A7*A3)/Z15*(Z17*Z7-Z1*Z16) 1+(A10*A1-A7*A4)/Z15*(Z1*Z13-Z17*Z8))/DEL C1=-A1*((A8*A1-A7*A2)/Z15*(Z13*Z4-Z1*Z16) 1+(A9*A1-A7*A3)/Z15*(Z6*Z16-Z17*Z4) 1+(A10*A1-A7*A4)/Z15*(Z17*Z1-Z6*Z13))/DEL GO TO 15 C**** EXPONENTIAL 14 TMP =A11*A2-A1*A1 E1 =(A11*ALXY-A1*ALL7)/TMP F1 =(ALL7*A2-A1*ALXY)/TMP GO TO 15 C**** POWER 16 TMP =A11*ALL2-ALL1*ALL1 E1 =(A11*ALLXY-ALL1*ALL7)/TMP F1 =(ALL7*ALL2-ALL1*ALLXY)/TMP 15 CALL OUT RETURN END SUBROUTINE GRAPH C INTEGER*2 IPERC(5),ICOUNT(6) LOGICAL PENPLT C C LOGICAL SCRN COMMON IBUF(500),Z(502) COMMON /A/ C1,E1,F1,D1,YE(502),IORDA,I2,Y(202),A(11),X(202) COMMON /B/ X1,X2,X3,X4,X5,X6,Y1,Y2,Y3,Y4,SCRN INTEGER TYPFIT(12,8) DATA TYPFIT /'LI','NE','AR',' T','HR','OU','GH',' O','RI','GI', * 'N ',' ', 1 'LI','NE','AR',' W','IT','H ','IN','TE','RC','EP', * 'T ',' ', 2 'QU','AD','RA','TI','C ','TH','RO','UG','H ','OR', * 'IG','IN', 3 'QU','AD','RA','TI','C ','WI','TH',' I','NT','ER', * 'CE','PT', 4 'CU','BI','C ','TH','RO','UG','H ','OR','IG','IN', * ' ',' ', 5 'CU','BI','C ','WI','TH',' I','NT','ER','CE','PT', * ' ',' ', 6 'EX','PO','NE','NT','IA','L ','(L','OG','-L','IN', * 'EA','R)', 7 'PO','WE','R ','(L','OG','-L','OG',') ',' ',' ', * ' ',' '/ DATA ICOUNT /'IN','DE','PE','ND','EN','T '/ DATA IPERC/'DE','PE','ND','EN','T '/ 6 CALL PLOTS(IBUF,500,1) CALL PLOT(1.5,1.5,-3) CALL SCALE(X(1),25.0,I2,1) XMIN=X(I2+1) XDEL=X(I2+2) CALL AXIS(0.,0.,ICOUNT,-11,25.0,0.0,X(I2+1),X(I2+2)) CALL SCALE(Y(1),19.0,I2,1) YMIN=Y(I2+1) YDEL=Y(I2+2) CALL AXIS(0.,0.,IPERC,10,19.0,90.,Y(I2+1),Y(I2+2)) CALL LINE(X,Y,I2,1,-1,4) C GET MAXIMUM OF X XMAX=.1E-36 DO 14 I=1,I2 IF(X(I).GT.XMAX)XMAX=X(I) 14 CONTINUE C COMPUTE MAXIMUM AXIS VALUE OF Y YMAX =YMIN+7.6*YDEL C DIF=(XMAX-XMIN)/500. Z(1)=XMAX YE(1)=F1+E1*XMAX+D1*XMAX**2+C1*XMAX**3 IF (IORDA.EQ.7) YE(1)=F1*E1**XMAX IF (IORDA.EQ.8) YE(1)=F1*XMAX**E1 ZZ=XMIN J=501 DO 15 I=1,499 J=J-1 ZZ=ZZ+DIF IF (ZZ.GT.XMAX) ZZ=XMAX 3 Z(J)=ZZ IF (IORDA.LT.7) GO TO 10 IF (IORDA.EQ.8) GO TO 8 C EXPONENTIAL YE(J) =F1*E1**ZZ GO TO 11 C POWER 8 YE(J) =F1*ZZ**E1 GO TO 11 C POLYNOMIAL 10 YE(J)=F1+E1*ZZ+D1*ZZ*ZZ+C1*ZZ**3 C C MAKE SURE Y VALUE IS WITHIN BOUNDS C 11 IF (YE(J).GE.YMIN.AND.YE(J).LE.YMAX) GO TO 15 C C Y VALUE IS OUT OF BOUNDS C IF (I.NE.1) GO TO 12 C C VERY FIRST VALUE. INCREMENT ZZ AND DO PARTIAL REITERATION C ZZ =ZZ+DIF GO TO 3 C C C MAKE X,Y VALUES SAME AS THAT OF PREVIOUS CELLS C 12 YE(J) =YE(J+1) Z(J) =Z(J+1) C 15 CONTINUE C C C NOW SEE IF THE PRESENT ENDPOINT VALUE IS WITHIN BOUNDS C IF (YE(1).GE.YMIN.AND.YE(1).LE.YMAX) GO TO 20 C OUT OF BOUND. SET IT TO ITS NEIGHBOUR'S VALUES YE(1) =YE(2) Z(1) =Z(2) 20 YE(501)=YMIN YE(502)=YDEL Z(501)=XMIN Z(502)=XDEL CALL LINE(Z,YE,500,1,0,0) C NOW PLOT THE TYPE-OF-FIT LABEL. CALL SYMBOL (5.,18.,0.75,TYPFIT(1,IORDA),0.,24) CALL PLOT (-1.5,-1.5,3) CALL PLOT (0.,0.,999) RETURN END SUBROUTINE OUT C INTEGER TYPE (12,8) C COMMON /A/C1,E1,F1,D1,YE(502),IORDA,I2,Y(202),A(11),X(202) C DATA TYPE /'LI','NE','AR',' T','HR','OU','GH',' O','RI','GI', * 'N ',' ', 1 'LI','NE','AR',' W','IT','H ','IN','TE','RC','EP', * 'T ',' ', 2 'QU','AD','RA','TI','C ','TH','RO','UG','H ','OR', * 'IG','IN', 3 'QU','AD','RA','TI','C ','WI','TH',' I','NT','ER', * 'CE','PT', 4 'CU','BI','C ','TH','RO','UG','H ','OR','IG','IN', * ' ',' ', 5 'CU','BI','C ','WI','TH',' I','NT','ER','CE','PT', * ' ',' ', 6 'EX','PO','NE','NT','IA','L ','(L','OG','-L','IN', * 'EA','R)', 7 'PO','WE','R ','(L','OG','-L','OG',') ',' ',' ', * ' ',' '/ C WRITE (2,200) (TYPE(I,IORDA),I=1,12) 200 FORMAT (20H0LEAST SQUARES FIT: ,12A2/1H ,43(1H-)) YBAR =A(7)/FLOAT(I2) GO TO (1,2,3,4,8,9,10,12),IORDA 1 WRITE(2,210)E1 210 FORMAT (5H0Y = ,E14.7,3H *X ) GO TO 5 2 WRITE(2,211)F1,E1 211 FORMAT (5H0 Y = ,E14.7,2H +,E14.7,3H *X ) GO TO 5 3 WRITE(2,212)E1,D1 212 FORMAT (5H0 Y = ,E14.7,5H *X +,E14.7,6H *X**2 ) GO TO 5 4 WRITE(2,213)F1,E1,D1 213 FORMAT (5H0Y = ,E14.7,2H +,E14.7,5H *X +,E14.7,6H *X**2 ) GO TO 5 8 WRITE(2,214)E1,D1,C1 214 FORMAT (5H0Y = ,E14.7,5H *X +,E14.7,8H *X**2 +,E14.7,6H +X**3 ) GO TO 5 9 WRITE(2,215)F1,E1,D1,C1 215 FORMAT (5H0Y = ,E14.7,2H +,E14.7,5H *X +,E14.7,8H *X**2 +,E14.7, 1 6H *X**3) GO TO 5 10 F1 =EXP(F1) WRITE (2,216) F1,E1 216 FORMAT (5H0Y = ,E14.7,6H *EXP(,E14.7,4H *X) ) E1 =EXP(E1) XK2 =0 XK3 =0 DO 11 I=1,I2 YE(I) =F1*E1**X(I) YEI =YE(I) YI =Y(I) XK1 =(YI-YEI)*(YI-YEI) XK2 =XK2+XK1 11 XK3 =XK3+(YI-YBAR)*(YI-YBAR) GO TO 14 12 F1 =EXP(F1) WRITE (2,217) F1,E1 217 FORMAT (5H0Y = ,E14.7,5H *X**,E14.7 ) XK2 =0 XK3 =0 DO 13 I=1,I2 YE(I) =F1*X(I)**E1 YEI =YE(I) YI =Y(I) XK1 =(YI-YEI)*(YI-YEI) XK2 =XK2+XK1 XK3 =XK3+(YI-YBAR)*(YI-YBAR) 13 CONTINUE GO TO 14 5 XK2=0 XK3 =0 DO 6 I=1,I2 YE(I)=C1*X(I)**3+D1*X(I)**2+E1*X(I)+F1 YI =Y(I) YEI =YE(I) XK1 =(YI-YEI)*(YI-YEI) XK3 =XK3+(YI-YBAR)*(YI-YBAR) 6 XK2=XK2+XK1 14 XI2=I2 S1=SQRT(XK2/XI2) S2=XI2*100.*S1/A(7) WRITE(2,228) 228 FORMAT (11H0 ACTUAL X,6X,8HACTUAL Y,4X,12HESTIMATED Y,2X, 1 12HPERCENT DEV/1H ,4(12(1H-),2X) ) DO 7 I=1,I2 IF (ABS(Y(I)).GT.1.E-5) GO TO 229 WRITE (2,220) X(I),Y(I),YE(I) 220 FORMAT (1H ,3(F10.3,4X),2X,8(1H*)) GO TO 7 229 DIX =(YE(I)-Y(I))/Y(I)*100 WRITE (2,202) X(I),Y(I),YE(I),DIX 7 CONTINUE 202 FORMAT (1H ,4(F10.3,4X)) WRITE(2,230)S1 230 FORMAT (23H0STANDARD ERROR OF EST.,2X,1H=,F6.2) WRITE(2,203)S2 203 FORMAT (26H RELATIVE STANDARD ERROR =,F5.1,1H%) CODET =1.-(XK2/XK3) WRITE (2,222) CODET 222 FORMAT (26H COEFF. OF DETERMINATION =,F6.2) RETURN END