Repatch (from linux) of OSX IRAF

1 files changed, 101 insertions, 0 deletions

diff --git a/math/ieee/chap1/fftohm.f b/math/ieee/chap1/fftohm.f
new file mode 100644
index 00000000..6e3ecc9d
--- /dev/null
+++ b/math/ieee/chap1/fftohm.f
@@ -0,0 +1,101 @@
+c
+c-----------------------------------------------------------------------
+c subroutine: fftohm
+c compute dft for real, n-point, odd harmonic sequences using an
+c n/2 point fft
+c odd harmonic means x(2*k)=0, all k where x(k) is the dft of x(m)
+c note: index m is sequence index--not fortran index
+c-----------------------------------------------------------------------
+c
+      subroutine fftohm(x, n)
+      dimension x(1)
+c
+c x = real array which on input contains the first n/2 points of the
+c     input
+c     on output x contains the n/4 complex values of the odd
+c     harmonics of the input--stored in the sequence re(x(1)),im(x(1)),
+c     re(x(2)),im(x(2)),...
+c ****note: x must be dimensioned to size n/2+2 for fft routine
+c n = true size of x sequence
+c
+c first compute real(x(1)) and real(x(n/2-1)) separately
+c also simultaneously multiply original sequence by sin(twopi*(m-1)/n)
+c sin and cos are computed recursively
+c
+c
+c for n = 2, assume x(1)=x0, x(2)=-x0, compute dft directly
+c
+      if (n.gt.2) go to 10
+      x(1) = 2.*x(1)
+      x(2) = 0.
+      return
+  10  twopi = 8.*atan(1.0)
+      tpn = twopi/float(n)
+c
+c compute x1=real(x(1)) and x2=imaginary(x(n/2-1))
+c x(n) = x(n)*4.*sin(twopi*(i-1)/n)
+c
+      t1 = 0.
+c
+c cosd and sind are multipliers for recursion for sin and cos
+c cosi and sini are initial conditions for recursion for sin and cos
+c
+      cosd = cos(tpn*2.)
+      sind = sin(tpn*2.)
+      cosi = 1.
+      sini = 0.
+      no2 = n/2
+      do 20 i=1,no2,2
+        t = x(i)*cosi
+        x(i) = x(i)*4.*sini
+        temp = cosi*cosd - sini*sind
+        sini = cosi*sind + sini*cosd
+        cosi = temp
+        t1 = t1 + t
+  20  continue
+c
+c reset initial conditions (cosi,sini) for new recursion
+c
+      cosi = cos(tpn)
+      sini = sin(tpn)
+      t2 = 0.
+      do 30 i=2,no2,2
+        t = x(i)*cosi
+        x(i) = x(i)*4.*sini
+        temp = cosi*cosd - sini*sind
+        sini = cosi*sind + sini*cosd
+        cosi = temp
+        t2 = t2 + t
+  30  continue
+      x1 = 2.*(t1+t2)
+      x2 = 2.*(t1-t2)
+c
+c take n/2 point (real) fft of preprocessed sequence x
+c
+      call fast(x, no2)
+c
+c for n = 4--skip recursion and initial conditions
+c
+      if (n.eq.4) go to 50
+c
+c initial conditions for recursion
+c
+      x(2) = -x(1)/2.
+      x(1) = x1
+c
+c for n = 8, skip recursion
+c
+      if (n.eq.8) go to 50
+c
+c unscramble y(k) using recursion formula
+c
+      nind = no2 - 2
+      do 40 i=3,nind,2
+        t = x(i)
+        x(i) = x(i-2) + x(i+1)
+        x(i+1) = x(i-1) - t
+  40  continue
+  50  x(no2) = x(no2+1)/2.
+      x(no2-1) = x2
+      return
+      end