Repatch (from linux) of OSX IRAF

1 files changed, 116 insertions, 0 deletions

diff --git a/math/ieee/chap1/fft842.f b/math/ieee/chap1/fft842.f
new file mode 100644
index 00000000..37a3a651
--- /dev/null
+++ b/math/ieee/chap1/fft842.f
@@ -0,0 +1,116 @@
+c
+c-----------------------------------------------------------------------
+c subroutine:  fft842
+c fast fourier transform for n=2**m
+c complex input
+c-----------------------------------------------------------------------
+c
+      subroutine fft842(in, n, x, y)
+c
+c this program replaces the vector z=x+iy by its  finite
+c discrete, complex fourier transform if in=0.  the inverse transform
+c is calculated for in=1.  it performs as many base
+c 8 iterations as possible and then finishes with a base 4 iteration
+c or a base 2 iteration if needed.
+c
+c the subroutine is called as subroutine fft842 (in,n,x,y).
+c the integer n (a power of 2), the n real location array x, and
+c the n real location array y must be supplied to the subroutine.
+c
+      dimension x(2), y(2), l(15)
+      common /con2/ pi2, p7
+      equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
+     *    (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
+     *    (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
+     *    (l1,l(15))
+c
+c
+c iw is a machine dependent write device number
+c
+      iw = i1mach(2)
+c
+      pi2 = 8.*atan(1.)
+      p7 = 1./sqrt(2.)
+      do 10 i=1,15
+        m = i
+        nt = 2**i
+        if (n.eq.nt) go to 20
+  10  continue
+      write (iw,9999)
+9999  format (35h n is not a power of two for fft842)
+      stop
+  20  n2pow = m
+      nthpo = n
+      fn = nthpo
+      if (in.eq.1) go to 40
+      do 30 i=1,nthpo
+        y(i) = -y(i)
+  30  continue
+  40  n8pow = n2pow/3
+      if (n8pow.eq.0) go to 60
+c
+c radix 8 passes,if any.
+c
+      do 50 ipass=1,n8pow
+        nxtlt = 2**(n2pow-3*ipass)
+        lengt = 8*nxtlt
+        call r8tx(nxtlt, nthpo, lengt, x(1), x(nxtlt+1), x(2*nxtlt+1),
+     *      x(3*nxtlt+1), x(4*nxtlt+1), x(5*nxtlt+1), x(6*nxtlt+1),
+     *      x(7*nxtlt+1), y(1), y(nxtlt+1), y(2*nxtlt+1), y(3*nxtlt+1),
+     *      y(4*nxtlt+1), y(5*nxtlt+1), y(6*nxtlt+1), y(7*nxtlt+1))
+  50  continue
+c
+c is there a four factor left
+c
+  60  if (n2pow-3*n8pow-1) 90, 70, 80
+c
+c go through the base 2 iteration
+c
+c
+  70  call r2tx(nthpo, x(1), x(2), y(1), y(2))
+      go to 90
+c
+c go through the base 4 iteration
+c
+  80  call r4tx(nthpo, x(1), x(2), x(3), x(4), y(1), y(2), y(3), y(4))
+c
+  90  do 110 j=1,15
+        l(j) = 1
+        if (j-n2pow) 100, 100, 110
+ 100    l(j) = 2**(n2pow+1-j)
+ 110  continue
+      ij = 1
+      do 130 j1=1,l1
+      do 130 j2=j1,l2,l1
+      do 130 j3=j2,l3,l2
+      do 130 j4=j3,l4,l3
+      do 130 j5=j4,l5,l4
+      do 130 j6=j5,l6,l5
+      do 130 j7=j6,l7,l6
+      do 130 j8=j7,l8,l7
+      do 130 j9=j8,l9,l8
+      do 130 j10=j9,l10,l9
+      do 130 j11=j10,l11,l10
+      do 130 j12=j11,l12,l11
+      do 130 j13=j12,l13,l12
+      do 130 j14=j13,l14,l13
+      do 130 ji=j14,l15,l14
+        if (ij-ji) 120, 130, 130
+ 120    r = x(ij)
+        x(ij) = x(ji)
+        x(ji) = r
+        fi = y(ij)
+        y(ij) = y(ji)
+        y(ji) = fi
+ 130    ij = ij + 1
+      if (in.eq.1) go to 150
+      do 140 i=1,nthpo
+        y(i) = -y(i)
+ 140  continue
+      go to 170
+ 150  do 160 i=1,nthpo
+        x(i) = x(i)/fn
+        y(i) = y(i)/fn
+ 160  continue
+ 170  return
+      end