Trailing-Edge
-
PDP-10 Archives
-
decus_20tap1_198111
-
decus/20-0025/orpol.for
There is 1 other file named orpol.for in the archive. Click here to see a list.
00010 C ORPOL.SRC OLS 148 LINES
00020 C
00030 COMMON Y(100),X(100),RHO(100),POLY(100),ALPHA(100),
00040 +BETA(30),COEF(30),P0(100),P1(100),P2(100)
00050 COMMON LL(30)
00060 DIMENSION IFMT(3)
00070 9901 FORMAT (A1)
00080 9902 FORMAT (3I)
00090 9910 FORMAT (10I)
00100 9900 FORMAT (10(F))
00110 TYPE 9899
00120 9899 FORMAT(1H ,'ORTHOGONAL POLYNOMIAL CURVE FITTING')
00130 XTEND=1.E37
00140 IYES='Y'
00150 1 TYPE 9898
00160 9898 FORMAT(' DO YOU WISH TO USE A DATA FILE, TYPE YES
00170 + OR NO? '$)
00180 ACCEPT 9901,IY
00190 IF (IYES-IY) 3, 150, 3
00200 3 IF (IY.EQ.'S') GO TO 160
00210 TYPE 9897
00220 9897 FORMAT(1H ,'TYPE NUMBER OF POINTS, MAXIMUM DEGREE '$)
00230 ACCEPT 9902,N,MAX
00240 MP1=MAX+1
00250 TYPE 9896
00260 9896 FORMAT(1H ,'TYPE IN DEPENDENT DATA'/1H )
00270 TYPE 1111,N,I
00280 1111 FORMAT(' N=',2I)
00290 CALL FREFLD(N,Y)
00300 TYPE 9890
00310 9890 FORMAT (1H ,'TYPE IN INDEPENDENT DATA'/1H )
00320 CALL FREFLD(N,X)
00330 TYPE 9895
00340 9895 FORMAT(1H ,'TYPE IN WEIGHTS'/1H )
00350 CALL FREFLD(N,RHO)
00360 4 FN=0.
00370 TA=0.
00380 BA=0.
00390 YSUM=0.
00400 DO 5 I=1,N
00410 FN=FN+RHO(I)
00420 5 YSUM=YSUM+Y(I)*RHO(I)
00430 YMEAN=YSUM/FN
00440 DO 10 I=1,N
00450 TA=TA+RHO(I)*X(I)
00460 10 BA=BA+RHO(I)
00470 ALPHA(1)=TA/BA
00480 BETA(I)=0.
00490 K=1
00500 SS=0.
00510 XPY=0.
00520 DO 20 I=1,N
00530 P0(I)=1.
00540 P1(I)=P0(I)*X(I)-ALPHA(1)*P0(I)
00550 SS=SS+P1(I)*P1(I)*RHO(I)
00560 Y(I)=Y(I)-YMEAN
00570 20 XPY=XPY+P1(I)*Y(I)*RHO(I)
00580 COEF(1)=XPY/SS
00590 SSR=COEF(1)*XPY
00600 TYPE 9894
00610 9894 FORMAT(1H ,'DEPENDENT DATA MEAN')
00620 TYPE 1001,YMEAN
00630 1001 FORMAT(1PE14.6)
00640 TYPE 9905
00650 9905 FORMAT (1H ,' DEGREE ALPHA BETA CO
00660 +EFF SSR')
00670 TYPE 1000,K,ALPHA(1),BETA(1),COEF(1),SSR
00680 1000 FORMAT(5X,I3,1P4E14.6)
00690 IF (MAX-1) 25, 61, 25
00700 25 DO 60 K=2,MAX
00710 BB=BA
00720 BA=0.
00730 TA=0.
00740 TB=0.
00750 DO 30 I=1,N
00760 TA=TA+RHO(I)*P1(I)*P1(I)*X(I)
00770 BA=BA+RHO(I)*P1(I)*P1(I)
00780 30 TB=TB+RHO(I)*X(I)*P1(I)*P0(I)
00790 ALPHA(K)=TA/BA
00800 BETA(K)=TB/BB
00810 SS=0.
00820 XPY=0.
00830 DO 40 I=1,N
00840 P2(I)=X(I)*P1(I)-ALPHA(K)*P1(I)-BETA(K)*P0(I)
00850 SS=SS+P2(I)*P2(I)*RHO(I)
00860 40 XPY=XPY+P2(I)*Y(I)*RHO(I)
00870 COEF(K)=XPY/SS
00880 SSR=COEF(K)*XPY
00890 TYPE 1000,K,ALPHA(K),BETA(K),COEF(K),SSR
00900 DO 50 J=1,N
00910 P0(J)=P1(J)
00920 50 P1(J)=P2(J)
00930 60 CONTINUE
00940 61 CONTINUE
00950 62 TYPE 9893
00960 9893 FORMAT(1H ,'DESIRED NUMBER OF POLYNOMIALS TO TRY? '$)
00970 ACCEPT 9902,M
00980 IF (M) 130,64,64
00990 64 TYPE 9892
01000 9892 FORMAT(1H ,'WHICH ONES? '$)
01010 ACCEPT 9910,(LL(I),I=1,M)
01020 65 TYPE 9891
01030 9891 FORMAT(1H ,'INPUT? '$)
01040 ACCEPT 9900,XT
01050 IF (XT-XTEND) 70, 62, 62
01060 70 POLY(1)=1.
01070 POLY(2)=XT-ALPHA(1)
01080 IF (MAX-1) 80, 95, 80
01090 80 DO 90 K=3,MP1
01100 90 POLY(K)=XT*POLY(K-1)-ALPHA(K-1)*POLY(K-1)-BETA(K-1)
01110 +*POLY(K-2)
01120 95 PRED=YMEAN
01130 DO 100 I=1,M
01140 J=LL(I)
01150 100 PRED=PRED+COEF(J)*POLY(J+1)
01160 TYPE 3009,PRED
01170 3009 FORMAT(16H PREDICTED VALUE ,1PE14.6)
01180 120 GO TO 65
01190 130 GO TO 1
01200 150 PRINT 9903
01210 9903 FORMAT (' FILE NAME? '$)
01220 ACCEPT 9904,FNM
01230 9904 FORMAT (A5)
01240 CALL IFILE (1,FNM)
01250 9990 IFMT(1)='('
01260 IFMT(3)='F)'
01270 READ (1,9902) N,MAX,NPL
01280 CALL ITOA(NPL,NEW)
01290 IFMT(2)=NEW
01300 NI=(N/NPL)+(MOD(N,NPL)/MOD(N,NPL))
01310 II=1
01320 DO 177 J=1,NI
01330 151 READ (1,IFMT)(Y(I),I=II,II+NPL-1)
01340 II=II+NPL
01350 177 CONTINUE
01360 II=1
01370 DO 155 J=1,NI
01380 152 READ (1,IFMT)(X(I),I=II,II+NPL-1)
01390 II=II+NPL
01400 155 CONTINUE
01410 II=1
01420 DO 159 J=1,NI
01430 153 READ (1,IFMT)(RHO(I),I=II,II+NPL-1)
01440 II =II+NPL
01450 159 CONTINUE
01460 MP1=MAX+1
01470 GO TO 4
01480 160 CALL EXIT
01490 END