Trailing-Edge
-
PDP-10 Archives
-
decuslib20-01
-
decus/20-0025/forir.for
There is 1 other file named forir.for in the archive. Click here to see a list.
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