Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-01 - 43,50144/orpol.src
There are no other files named orpol.src in the archive.
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