File: fourm.x

package info (click to toggle)
iraf-rvsao 2.8.3-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 16,456 kB
  • sloc: ansic: 963; lisp: 651; fortran: 397; makefile: 27
file content (85 lines) | stat: -rw-r--r-- 1,849 bytes parent folder | download
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
# File rvsao/Util/fourm.x
# September 16, 1991
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# From Press, et al., Numerical Recipes FOUR1
# Modified by Dan Fabricant, Center for Astrophysics

#--- Perform forward (isign=1) or reverse (isign=-1) Fourier transform on vector

procedure fourm (xdata,nin,isign)

real	xdata[ARB]	# Complex vector to transform (result returned)
int	nin		# Log base 2 of length of vector
int	isign		# 1=forward, -1=reverse

double	wr,wi,wpr,wpi,wtemp,theta
real	xr, xi, tempi, tempr
int	i,j,m,mmax,n,nn,istep

begin
	nn = 2 ** nin
	n = 2 * nn
	j = 1

#  Exchange the two complex numbers
	do i = 1 ,n, 2 {
	    if (j > i) {
		tempr = xdata[j]
		tempi = xdata[j+1]
		xdata[j] = xdata[i]
		xdata[j+1] = xdata[i+1]
		xdata[i] = tempr
		xdata[i+1] = tempi
		}
	    m = n / 2

	    while (m >= 2 && j > m) {
		j = j - m
		m = m / 2
		}
	    j = j+m
	    }

#  Danielson-Lanczos section executed nin times
	mmax = 2
	while (n > mmax) {
	    istep = 2 * mmax
	    theta = 6.28318530717959d0 / (isign * mmax)
	    wpr = -2.d0 * ((sin (0.5d0 * theta)) ** 2)
	    wpi = sin (theta)
	    wr = 1.d0
	    wi = 0.d0

	    do m = 1, mmax, 2 {
		do i = m, n, istep {
		    j = i + mmax

		#  Danielson-Lanczos formula
		    xr = real (wr)
		    xi = real (wi)
		    tempr = (xr * xdata[j]) - (xi * xdata[j+1])
		    tempi = (xr * xdata[j+1]) + (xi * xdata[j])
		    xdata[j] = xdata[i] - tempr
		    xdata[j+1] = xdata[i+1] - tempi

		    xdata[i] = xdata[i] + tempr
		    xdata[i+1] = xdata[i+1] + tempi
		    }
	    # Trigonometric recurrence
		wtemp = wr
		wr = (wr * wpr) - (wi * wpi) + wr
		wi = (wi * wpr) + (wtemp * wpi) + wi
		}
	    mmax = istep
	    }

# Added to make plug compatible with previously used subroutine (FFT2C)

	if (isign < 0) {
	    do i = 1, n {
		xdata[i] = xdata[i] / nn
		}
	    }

	return
end