Repatch (from linux) of OSX IRAF

1 files changed, 488 insertions, 0 deletions

diff --git a/math/ieee/chap1/radix4.f b/math/ieee/chap1/radix4.f
new file mode 100644
index 00000000..cd703305
--- /dev/null
+++ b/math/ieee/chap1/radix4.f
@@ -0,0 +1,488 @@
+c
+c-----------------------------------------------------------------------
+c subroutine:  radix4
+c computes forward or inverse complex dft via radix-4 fft.
+c uses autogen technique to yield time efficient program.
+c-----------------------------------------------------------------------
+c
+        subroutine radix4(mm,iflag,jflag)
+c
+c       input:
+c             mm = power of 4 (i.e., n = 4**mm complex point transform)
+c             (mm.ge.2 and mm.le.5)
+c
+c          iflag = 1 on first pass for given n
+c                = 0 on subsequent passes for given n
+c
+c          jflag = -1 for forward transform
+c                = +1 for inverse transform
+c
+c       input/output:
+c              a = array of dimensions 2*n with real and imaginary parts
+c                  of dft input/output in odd, even array components.
+c
+c       for optimal time efficiency, common is used to pass arrays.
+c       this means that dimensions of arrays a, ix, and t can be
+c       modified to reflect maximum value of n = 4**mm to be used. note
+c       that array "ix" is also dimensioned in subroutine "rad4sb".
+c
+c    i.e.,    a(    ) ix(  ) t(   )
+c
+c       m =2      32     38     27
+c       m<=3     128    144    135
+c       m<=4     512    658    567
+c       m<=5    2048   2996   2295
+c
+       dimension a(2048),ix(2996),t(2295)
+       dimension nfac(11),np(209)
+       common ntypl,kkp,index,ixc
+       common /aa/a
+       common /xx/ix
+c
+c      check for mm<2 or mm>5
+c
+       if(mm.lt.2.or.mm.gt.5)stop
+c
+c      initialize on first pass """"""""""""""""""""""""""""""""""""""""
+c
+       if(iflag.eq.1) go to 9999
+c
+c      fast fourier transform start ####################################
+c
+8885   kspan=2*4**mm
+       if(jflag.eq.1) go to 8887
+c
+c      conjugate data for forward transform
+c
+       do 8886 j=2,n2,2
+8886   a(j)=-a(j)
+       go to 8889
+c
+c      multiply data by n**(-1) if inverse transform
+c
+8887   do 8888 j=1,n2,2
+       a(j)=a(j)*xp
+8888   a(j+1)=a(j+1)*xp
+8889   i=3
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8),it
+c***********************************************************************
+c
+c                                    8 multiply butterfly
+c
+1      kk=ix(i)
+c
+11     k1=kk+kspan
+       k2=k1+kspan
+       k3=k2+kspan
+c
+       akp=a(kk)+a(k2)
+       akm=a(kk)-a(k2)
+       ajp=a(k1)+a(k3)
+       ajm=a(k1)-a(k3)
+       a(kk)=akp+ajp
+c
+       bkp=a(kk+1)+a(k2+1)
+       bkm=a(kk+1)-a(k2+1)
+       bjp=a(k1+1)+a(k3+1)
+       bjm=a(k1+1)-a(k3+1)
+       a(kk+1)=bkp+bjp
+c
+       bjp=bkp-bjp
+c
+       a(k2+1)=(akp+bjp-ajp)*c707
+       a(k2)=a(k2+1)+bjp*cm141
+c
+       bkp=bkm+ajm
+       akp=akm-bjm
+c
+       ac0=(akp+bkp)*c924
+       a(k1+1)=ac0+akp*cm541
+       a(k1)  =ac0+bkp*cm131
+c
+       bkm=bkm-ajm
+       akm=akm+bjm
+c
+       ac0=(akm+bkm)*c383
+       a(k3+1)=ac0+akm*c541
+       a(k3)  =ac0+bkm*cm131
+c
+       i=i+1
+       kk=ix(i)
+       if (kk) 111,111,11
+111    i=i+2
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8), it
+c***********************************************************************
+c
+c                                    4 multiply butterfly
+c
+2      kk=ix(i)
+c
+22     k1=kk+kspan
+       k2=k1+kspan
+       k3=k2+kspan
+c
+       akp=a(kk)+a(k2)
+       akm=a(kk)-a(k2)
+       ajp=a(k1)+a(k3)
+       ajm=a(k1)-a(k3)
+       a(kk)=akp+ajp
+c
+       bkp=a(kk+1)+a(k2+1)
+       bkm=a(kk+1)-a(k2+1)
+       bjp=a(k1+1)+a(k3+1)
+       bjm=a(k1+1)-a(k3+1)
+       a(kk+1)=bkp+bjp
+       a(k2)=-bkp+bjp
+       a(k2+1)=akp-ajp
+c
+       bkp=bkm+ajm
+c
+       a(k1+1)=(bkp+akm-bjm)*c707
+       a(k1)=a(k1+1)+bkp*cm141
+c
+       akm=akm+bjm
+c
+       a(k3+1)=(akm+ajm-bkm)*c707
+       a(k3)=a(k3+1)+akm*cm141
+c
+       i=i+1
+       kk=ix(i)
+       if (kk) 222,222,22
+222    i=i+2
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8), it
+c***********************************************************************
+c
+c                                    8 multiply butterfly
+c
+3      kk=ix(i)
+c
+33     k1=kk+kspan
+       k2=k1+kspan
+       k3=k2+kspan
+c
+       akp=a(kk)+a(k2)
+       akm=a(kk)-a(k2)
+       ajp=a(k1)+a(k3)
+       ajm=a(k1)-a(k3)
+       a(kk)=akp+ajp
+c
+       bkp=a(kk+1)+a(k2+1)
+       bkm=a(kk+1)-a(k2+1)
+       bjp=a(k1+1)+a(k3+1)
+       bjm=a(k1+1)-a(k3+1)
+       a(kk+1)=bkp+bjp
+c
+       ajp=akp-ajp
+c
+       a(k2+1)=(ajp+bjp-bkp)*c707
+       a(k2)=a(k2+1)+ajp*cm141
+c
+       bkp=bkm+ajm
+       akp=akm-bjm
+c
+       ac0=(akp+bkp)*c383
+       a(k1+1)=ac0+akp*c541
+       a(k1)  =ac0+bkp*cm131
+c
+       bkm=bkm-ajm
+       akm=akm+bjm
+c
+       ac0=(akm+bkm)*cm924
+       a(k3+1)=ac0+akm*c541
+       a(k3)  =ac0+bkm*c131
+c
+       i=i+1
+       kk=ix(i)
+       if (kk) 333,333,33
+333    i=i+2
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8), it
+c***********************************************************************
+c
+c                                    general 9 multiply butterfly
+c
+4      kk=ix(i)
+c
+44     k1=kk+kspan
+       k2=k1+kspan
+       k3=k2+kspan
+c
+       akp=a(kk)+a(k2)
+       akm=a(kk)-a(k2)
+       ajp=a(k1)+a(k3)
+       ajm=a(k1)-a(k3)
+       a(kk)=akp+ajp
+c
+       bkp=a(kk+1)+a(k2+1)
+       bkm=a(kk+1)-a(k2+1)
+       bjp=a(k1+1)+a(k3+1)
+       bjm=a(k1+1)-a(k3+1)
+       a(kk+1)=bkp+bjp
+c
+       ajp=akp-ajp
+       bjp=bkp-bjp
+c
+       j=ix(i+1)
+c
+       ac0=(ajp+bjp)*t(j+8)
+       a(k2+1)=ac0+ajp*t(j+6)
+       a(k2)  =ac0+bjp*t(j+7)
+c
+       bkp=bkm+ajm
+       akp=akm-bjm
+c
+       ac0=(akp+bkp)*t(j+5)
+       a(k1+1)=ac0+akp*t(j+3)
+       a(k1)  =ac0+bkp*t(j+4)
+c
+       bkm=bkm-ajm
+       akm=akm+bjm
+c
+       ac0=(akm+bkm)*t(j+2)
+       a(k3+1)=ac0+akm*t(j)
+       a(k3)  =ac0+bkm*t(j+1)
+c
+       i=i+2
+       kk=ix(i)
+       if (kk) 444,444,44
+444    i=i+2
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8), it
+c***********************************************************************
+c
+c                                    0 multiply butterfly
+c
+5      kk=ix(i)
+c
+55     k1=kk+kspan
+       k2=k1+kspan
+       k3=k2+kspan
+c
+       akp=a(kk)+a(k2)
+       akm=a(kk)-a(k2)
+       ajp=a(k1)+a(k3)
+       ajm=a(k1)-a(k3)
+       a(kk)=akp+ajp
+       a(k2)=akp-ajp
+c
+       bkp=a(kk+1)+a(k2+1)
+       bkm=a(kk+1)-a(k2+1)
+       bjp=a(k1+1)+a(k3+1)
+       bjm=a(k1+1)-a(k3+1)
+       a(kk+1)=bkp+bjp
+       a(k2+1)=bkp-bjp
+c
+       a(k3+1)=bkm-ajm
+       a(k1+1)=bkm+ajm
+       a(k3)=akm+bjm
+       a(k1)=akm-bjm
+c
+       i=i+1
+       kk=ix(i)
+       if (kk) 555,555,55
+555    i=i+2
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8), it
+c***********************************************************************
+c
+c                                       offset reduced
+c
+6      kspan=kspan/4
+       i=i+2
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8), it
+c***********************************************************************
+c
+c                                       bit reversal (shuffling)
+c
+7      ip1=ix(i)
+77     ip2=ix(i+1)
+       t1=a(ip2)
+       a(ip2)=a(ip1)
+       a(ip1)=t1
+       t1=a(ip2+1)
+       a(ip2+1)=a(ip1+1)
+       a(ip1+1)=t1
+       i=i+2
+       ip1=ix(i)
+       if (ip1) 777,777,77
+777    i=i+2
+       it=ix(i-1)
+       go to (1,2,3,4,5,6,7,8), it
+c***********************************************************************
+8      if(jflag.eq.1) go to 888
+c
+c      conjugate output if forward transform
+c
+       do 88 j=2,n2,2
+88     a(j)=-a(j)
+888    return
+c
+c      fast fourier transform ends #####################################
+c
+c      initialization phase starts. done only once
+c
+9999   ixc=1
+       n=4**mm
+       xp=n
+       xp=1./xp
+       ntot=n
+       n2=n*2
+       nspan=n
+       n1test=n/16
+       n2test=n/8
+       n3test=(3*n)/16
+       nspan4=nspan/4
+       ibase=0
+       isn=1
+       inc=isn
+       rad=8.0*atan(1.0)
+       pi=4.*atan(1.0)
+       c707=sin(pi/4.)
+       cm141=-2.*c707
+       c383=sin(pi/8.)
+       c924=cos(pi/8.)
+       cm924=-c924
+       c541=c924-c383
+       cm541=-c541
+       c131=c924+c383
+       cm131=-c131
+10     nt=inc*ntot
+       ks=inc*nspan
+       kspan=ks
+       jc=ks/n
+       radf=rad*float(jc)*.5
+       i=0
+c
+c      determine the factors of n
+c      all factors must be 4 for this version
+c
+       m=0
+       k=n
+15     m=m+1
+       nfac(m)=4
+       k=k/4
+20     if(k-(k/4)*4.eq.0) go to 15
+       kt=1
+       if(n.ge.256) kt=2
+       kspan0=kspan
+       ntypl=0
+c
+100    ndelta=kspan0/kspan
+       index=0
+       sd=radf/float(kspan)
+       cd=2.0*sin(sd)**2
+       sd=sin(sd+sd)
+       kk=1
+       i=i+1
+c
+c      transform for a factor of 4
+c
+       kspan=kspan/4
+       ix(ixc)=0
+       ix(ixc+1)=6
+       ixc=ixc+2
+c
+410    c1=1.0
+       s1=0.0
+420    k1=kk+kspan
+       k2=k1+kspan
+       k3=k2+kspan
+       if(s1.eq.0.0) go to 460
+430    if(kspan.ne.nspan4) go to 431
+       t(ibase+5)=-(s1+c1)
+       t(ibase+6)=c1
+       t(ibase+4)=s1-c1
+       t(ibase+8)=-(s2+c2)
+       t(ibase+9)=c2
+       t(ibase+7)=s2-c2
+       t(ibase+2)=-(s3+c3)
+       t(ibase+3)=c3
+       t(ibase+1)=s3-c3
+       ibase=ibase+9
+c
+431    kkp=(kk-1)*2
+       if(index.ne.n1test) go to 150
+       call rad4sb(1)
+       go to 5035
+150    if(index.ne.n2test) go to 160
+       call rad4sb(2)
+       go to 5035
+160    if(index.ne.n3test) go to 170
+       call rad4sb(3)
+       go to 5035
+170    call rad4sb(4)
+5035   kk=k3+kspan
+       if(kk.le.nt) go to 420
+440    index=index+ndelta
+       c2=c1-(cd*c1+sd*s1)
+       s1=(sd*c1-cd*s1)+s1
+       c1=c2
+       c2=c1*c1-s1*s1
+       s2=c1*s1+c1*s1
+       c3=c2*c1-s2*s1
+       s3=c2*s1+s2*c1
+       kk=kk-nt+jc
+       if(kk.le.kspan) go to 420
+       kk=kk-kspan+inc
+       if(kk.le.jc) go to 410
+       if(kspan.eq.jc) go to 800
+       go to 100
+460    kkp=(kk-1)*2
+       call rad4sb(5)
+5050   kk=k3+kspan
+       if(kk.le.nt) go to 420
+       go to 440
+c
+800    ix(ixc)=0
+       ix(ixc+1)=7
+       ixc=ixc+2
+c
+c      compute parameters to permute the results to normal order
+c      done in two steps
+c      permutation for square factors of n
+c
+       np(1)=ks
+       k=kt+kt+1
+       if(m.lt.k) k=k-1
+       j=1
+       np(k+1)=jc
+810    np(j+1)=np(j)/nfac(j)
+       np(k)=np(k+1)*nfac(j)
+       j=j+1
+       k=k-1
+       if(j.lt.k) go to 810
+       k3=np(k+1)
+       kspan=np(2)
+       kk=jc+1
+       k2=kspan+1
+       j=1
+c
+c      permutation for single variate transform
+c
+820    kkp=(kk-1)*2
+       k2p=(k2-1)*2
+       ix(ixc)=kkp+1
+       ix(ixc+1)=k2p+1
+       ixc=ixc+2
+       kk=kk+inc
+       k2=kspan+k2
+       if(k2.lt.ks) go to 820
+830    k2=k2-np(j)
+       j=j+1
+       k2=np(j+1)+k2
+       if(k2.gt.np(j)) go to 830
+       j=1
+840    if(kk.lt.k2) go to 820
+       kk=kk+inc
+       k2=kspan+k2
+       if(k2.lt.ks) go to 840
+       if(kk.lt.ks) go to 830
+       jc=k3
+       ix(ixc)=0
+       ix(ixc+1)=8
+       go to 8885
+       end