annotate libcruft/fftpack/cffti1.f @ 5103:e2ed74b9bfa0 after-gnuplot-split

[project @ 2004-12-28 02:43:01 by jwe]
author jwe
date Tue, 28 Dec 2004 02:43:01 +0000
parents 44ed237bdc1e
children 82be108cc558
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1644
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
1 subroutine cffti1 (n,wa,ifac)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
2 implicit double precision (a-h,o-z)
1645
44ed237bdc1e [project @ 1995-12-14 08:32:49 by jwe]
jwe
parents: 1644
diff changeset
3 dimension wa(*) ,ifac(*) ,ntryh(4)
1644
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
4 data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
5 nl = n
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
6 nf = 0
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
7 j = 0
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
8 101 j = j+1
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
9 if (j-4) 102,102,103
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
10 102 ntry = ntryh(j)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
11 go to 104
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
12 103 ntry = ntry+2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
13 104 nq = nl/ntry
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
14 nr = nl-ntry*nq
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
15 if (nr) 101,105,101
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
16 105 nf = nf+1
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
17 ifac(nf+2) = ntry
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
18 nl = nq
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
19 if (ntry .ne. 2) go to 107
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
20 if (nf .eq. 1) go to 107
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
21 do 106 i=2,nf
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
22 ib = nf-i+2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
23 ifac(ib+2) = ifac(ib+1)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
24 106 continue
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
25 ifac(3) = 2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
26 107 if (nl .ne. 1) go to 104
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
27 ifac(1) = n
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
28 ifac(2) = nf
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
29 tpi = 6.28318530717959d0
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
30 argh = tpi/dble(n)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
31 i = 2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
32 l1 = 1
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
33 do 110 k1=1,nf
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
34 ip = ifac(k1+2)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
35 ld = 0
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
36 l2 = l1*ip
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
37 ido = n/l2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
38 idot = ido+ido+2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
39 ipm = ip-1
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
40 do 109 j=1,ipm
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
41 i1 = i
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
42 wa(i-1) = 1.
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
43 wa(i) = 0.
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
44 ld = ld+l1
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
45 fi = 0.
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
46 argld = dble(ld)*argh
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
47 do 108 ii=4,idot,2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
48 i = i+2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
49 fi = fi+1.
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
50 arg = fi*argld
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
51 wa(i-1) = cos(arg)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
52 wa(i) = sin(arg)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
53 108 continue
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
54 if (ip .le. 5) go to 109
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
55 wa(i1-1) = wa(i-1)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
56 wa(i1) = wa(i)
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
57 109 continue
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
58 l1 = l2
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
59 110 continue
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
60 return
395bb6d6c096 [project @ 1995-12-14 08:32:12 by jwe]
jwe
parents:
diff changeset
61 end