forked from geodynamics/hc
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rick_fft.f90
196 lines (177 loc) · 4.92 KB
/
rick_fft.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
!
! fourier transform routines as used by rick_sh routines
! based on Rick O'Connell's subroutines, which are modified
! Numerical Recipes
!
! $Id: rick_fft.f90,v 1.4 2006/01/22 02:11:12 becker Exp becker $
!
subroutine rick_f90_cs2ab(rdata,n)
!
! Transforms spectral coefficients from cos-sin series to
! complex discrete fourier series. Function is real, and
! transformed by realft(rdata,n/2,1). Number of data points
! is n. Does not recover real component for frequency n/2.
use rick_module
implicit none
integer :: n,i
real(kind=cp), intent(inout),dimension(n) :: rdata
real(kind=cp) :: en
en = dfloat(n)
rdata(1) = rdata(1) * en
do i=3,n
rdata(i) = rdata(i)*en/2.0_cp
enddo
return
end subroutine rick_f90_cs2ab
subroutine rick_f90_ab2cs(rdata,n)
! Changes coefficients of complex spectrum of a real function
! transformed by realft.f to real coefficients of a series
! of C*cos(m*x)+S*sin(mx). Coefficients are ordered as
! C(0),S(0),C(1),S(1),C(2),...,C(n/2-1),S(n/2-1). This loses
! the real part of spectrum for frequency n/2.
! The number of data points is n, The call to realft is
! call realft(rdata,n/2,1)
!
use rick_module
implicit none
integer n,i
real(kind=cp), intent(inout), dimension(n) :: rdata
real(kind=cp) :: en
en=dfloat(n)
rdata(1)=rdata(1)/en
rdata(2)= 0.0_cp
do i=3,n
rdata(i)=rdata(i)/en*2.0_cp
enddo
return
end subroutine rick_f90_ab2cs
subroutine rick_f90_realft(rdata,n,isign)
!
! Calculates the fourier transform of 2*N real data points.
! Replaces data with the positive frequency half of the
! complex fourier transform. The real parts of the first
! and last frequency components are returned in data(1)
! and data(2) (i.e. for frequencies of zero and N/2). The
! other spectral components are given as complex pairs
! in data(3),data(4) etc. The inverse transform is obtained
! with ISIGN=-1, and dividing the data or result by N.
! Calls routine four1(data,n,isign) for FFT.
!
use rick_module
implicit none
integer, intent(in) :: n,isign
real(kind=cp), intent(inout), dimension(2*n+2) :: rdata
! local
real(kind=cp) :: c1,c2,h1r,h1i,h2r,h2i
double precision :: theta,wi,wpi,wpr,wr,wtemp
integer :: i,n2p3,ilim,i1,i2,i3,i4,n2
theta = pi / dble(n) !this is different from Numrec
wr = 1.0d0
wi = 0.0d0
c1 = 0.5_cp
n2 = 2*n
n2p3 = n2 + 3
if (isign.eq.1) then
c2=-0.5_cp
call rick_f90_four1(rdata,n,+1)
rdata(n2+1)=rdata(1)
rdata(n2+2)=rdata(2)
else
c2=0.5_cp
theta=-theta
rdata(n2+1)=rdata(2)
rdata(n2+2)=0.0_cp
rdata(2)=0.0_cp
endif
wpr=-2.0d0 * sin(0.5d0 * theta)**2
wpi=sin(theta)
ilim = n/2+1
do i=1,ilim
i1=2*i-1 ! 1,3,5,...
i2=i1+1
i3=n2p3-i2
i4=i3+1
h1r=c1*(rdata(i1)+rdata(i3))
h1i=c1*(rdata(i2)-rdata(i4))
h2r=-c2*(rdata(i2)+rdata(i4))
h2i= c2*(rdata(i1)-rdata(i3))
rdata(i1)= h1r+wr*h2r-wi*h2i
rdata(i2)= h1i+wr*h2i+wi*h2r
rdata(i3)= h1r-wr*h2r+wi*h2i
rdata(i4)=-h1i+wr*h2i+wi*h2r
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
enddo
if (isign.eq.1) then
rdata(2)=rdata(n2+1)
else
call rick_f90_four1(rdata,n,-1)
endif
return
end subroutine rick_f90_realft
subroutine rick_f90_four1(rdata,nn,isign)
!
! FFT routine from Numerical Recipes. Replaces data by
! its discrete fourier transform if isign=1, or by
! NN times its inverse transform if isign=-1. Array
! data is made up of NN complex numbers (2*NN pairs)
! and NN must be a power of 2. Spectral components
! are complex, and ordered from frequency zero to
! +-NN/2 to -1 in the standard fashion.
!
use rick_module
implicit none
integer, intent(in) :: nn,isign
real(kind=cp), intent(inout), dimension(2*nn+2) :: rdata
! local
real(kind=cp) :: tempr,tempi
! this should be double precision locally, regardless
double precision :: wr,wi,wpr,wpi,wtemp,theta
integer :: n,m,i,j,mmax,istep
n=2*nn
j=1
do i=1,n,2
if(j.gt.i)then
tempr=rdata(j)
tempi=rdata(j+1)
rdata(j)=rdata(i)
rdata(j+1)=rdata(i+1)
rdata(i)=tempr
rdata(i+1)=tempi
endif
m=n/2
1 if ((m.ge.2).and.(j.gt.m)) then
j=j-m
m=m/2
go to 1
endif
j=j+m
enddo
mmax=2
2 if (n.gt.mmax) then
istep=2*mmax
theta = 6.28318530717959d0/(isign*mmax)
wpr=-2.0d0 * sin(0.5d0 * theta)**2
wpi=sin(theta)
wr = 1.0d0
wi = 0.0d0
do m=1,mmax,2
do i=m,n,istep
j=i+mmax
tempr=wr*rdata(j)-wi*rdata(j+1)
tempi=wr*rdata(j+1)+wi*rdata(j)
rdata(j)=rdata(i)-tempr
rdata(j+1)=rdata(i+1)-tempi
rdata(i)=rdata(i)+tempr
rdata(i+1)=rdata(i+1)+tempi
enddo
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
enddo
mmax=istep
go to 2
endif
return
end subroutine rick_f90_four1