Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-01 - 43,50144/forir.f4
There are no other files named forir.f4 in the archive.
00010	        COMMON F(100),V(100),SN(100),CS(100),U(100)
00020	     1,A(100),B(100),SS(100),LL(100)
00021		TYPE 899
00022		TYPE 7001
00024	 7001	FORMAT(/' NUMBER OF DATA POINTS CANNOT
00025	     + EXCEED 100.'/' TO FINISH PREDICTED VALUES, 
00027	     + ENTER A'/' VALUE LARGER THAN 1.E10.')
00040	       ZERO=0.0
00060	       PI=3.14159265
00062	    1  TYPE 1111
00064	 1111  FORMAT(/'  TOTAL NUMBER OF DATA POINTS? ',$)
00066	      ACCEPT 56,LIM
00072	      IF(LIM)999,999,31
00074	   31 NROW = LIM/10
00076		LAST = LIM - 10*NROW
00078		TYPE 1031,NROW,LAST
00079		IF(NROW)34,34,32
00080	   32	L = 0
00082		DO 33 I = 1,NROW
00084		K = L + 1
00086		L = L + 10
00088		TYPE 1032,I
00090		ACCEPT 1033,(F(J),J=K,L)
00092	   33	CONTINUE
00094	   34	K = L + 1
00096		L = L + 10
00098		TYPE 1034
00100		ACCEPT 1033,(F(J),J=K,L)
00102		TYPE 1035
00104	 1031	FORMAT(/'  ENTER DATA ON',I3,' ROWS OF 10
00108	     +  POINTS EACH.'/'  ENTER REMAINING',I3,' POINTS
00110	     + ON LAST LINE.'/)
00112	 1032	FORMAT('  ROW',I3,' ? ',$)
00114	 1033	FORMAT(10F)
00116	 1034	FORMAT('  LAST   ? ',$)
00118	 1035	FORMAT('  THANK YOU'/)
00190	 60    N=LIM/2
00200	       FN=N
00210	       M=N-1
00220	       FLIM=(2.*FN-1.)/FN
00230	C BASIC SINES AND COSINES
00240	       V(1)=0.0
00250	       V(2)=1.0
00260	       CF=COS(PI/FN)
00270	       SF=SIN(PI/FN)
00280	       MP1=M+1
00290	       DO 80 K=3,MP1
00300	       V(K)=2.*CF*V(K-1)-V(K-2)
00310	       SN(K-1)=SF*V(K)
00320	       CS(K-1)=CF*V(K)-V(K-1)
00330	 80    CONTINUE
00340	       CS(1)=CF
00350	       SN(1)=SF
00360	 899   FORMAT(24H   FINITE FOURIER SERIES)
00370	C COEFFICIENT LOOP
00380	       SSUM=0.
00390	       DO 82 I=1,LIM
00400	 82    SSUM=F(I)+SSUM
00410	       AZERO=SSUM/FN
00420	       U(1)=0.
00430	       U(2)=F(LIM)
00440	       DO 95 K=1,M
00450	       DO 85 I=3,LIM
00460	       INDEX=LIM-I+2
00470	 85    U(I)=U(I-1)*2.*CS(K)-U(I-2)+F(INDEX)
00480	       A(K)=(CS(K)*U(LIM)-U(LIM-1)+F(1))/FN
00490	       B(K)=SN(K)*U(LIM)/FN
00500	 95    CONTINUE
00510	       SUMS=0.
00520	       DO 101 I=1,LIM
00530	 101   SUMS=SUMS+F(I)*F(I)
00540	C CALCULATE ERROR SUM OF SQUARES
00550	 107   TEMP=FN*AZERO*AZERO/2.
00560	       SSO=TEMP
00570	       SS(1)=(A(1)**2+B(1)**2)*FN
00580	       DO 108 K=2,M
00590	 108   SS(K)=(A(K)**2+B(K)**2)*FN
00600	C PRINT COEFFICIENTS AND ERROR SUM OF SQUARES
00610		TYPE 1001
00620	 1001  FORMAT(/'  HARMONIC   COS TERMS    SINE TERMS
00630	     1   ERROR SS REMOVED')
00650	       IOO=0
00660	       TYPE  99, IOO,AZERO,ZERO,SSO
00670	 99    FORMAT(8X,I2,1P3E14.6)
00680	       DO 109 I=1,M
00690	 109   TYPE  99, I,A(I),B(I),SS(I)
00700	111	TYPE 55
00701	55	FORMAT(/'  DESIRED NO. OF HARMONICS TO TRY? '$)
00720		ACCEPT 56,NH
00721	56	FORMAT (I )
00722		IF(NH-1-M)110,110,111
00730	110     IF(NH)130,130,112
00740	112	TYPE 57
00741	57	FORMAT (' WHICH ONES? '$)
00750	 115   ACCEPT 58,  (LL(I),I=1,NH)
00751	58	FORMAT (10I )
00760	118	TYPE 59
00761	59	FORMAT(/' INPUT? '$)
00770		ACCEPT 41,XT
00771	41	FORMAT(E)
00780	       IF (XT-1.E10) 120,111,111
00790	 120   TI=0.
00800	       DO 125 K=1,NH
00810	       J=LL(K)
00820	       IF (J) 123,122,123
00830	 122   TI=TI+AZERO/2.
00840	       GO TO 125
00850	 123   FJ=J
00855	      ANGLE = PI * FJ * XT / FN
00860	      TI = TI + A(J)*COS(ANGLE) + B(J)*SIN(ANGLE)
00880	 125   CONTINUE
00890	       TYPE  1025,TI
00900	 1025  FORMAT(16H+PREDICTED VALUE,1PE14.6)
00910	 	GO TO 118
00920	 130	GO TO 1
00930	999   CALL EXIT
00940	      END