......subroutine fft2(ar,ai,n) <br> 
c fft routine with 2 real input arrays rather than complex <br> 
c fft2 subroutine mults by exp(i\omega t) <br> 
c OUTPUT: f=0 in ar(1), f=f_N in ar(n/2+1) <br> 
c <br> 
c n is a power of two. <br> 
c For a real-valued time series, the FT of negative frequencies resides in the <br> 
c output array beyond the Nyquist f_N with the usual conjugation: <br> 
c <br> 
c ar(i)=ar(n+2-i), ai(i)=-ai(n+2-i) <br> 
c <br> 
c fft2 is NOT a unitary transform, mults the series by sqrt(n) <br> 
c the inverse FT can be effected by running fft2 on the conjugate of <br> 
c the FFT-expansion, then taking the the conjugate of the output, <br> 
c and dividing thru by N. to wit: <br> 
c <br> 
c assume Xr, Xi is in freq domain, xr, xi in the time domain <br> 
c <br> 
c (Xr,Xi)=fft2(xr,xi,N) <br> 
c (xr,-xi)=fft2(Xr,-Xi,N)/N <br> 
c <br> 
......implicit real*8 (a-h,o-z) <br> 
......implicit integer*4 (i-n) <br> 
......dimension ar(1),ai(1) <br> 
......mex=dlog(dble(float(n)))/.693147d0 <br> 
......nv2=n/2 <br> 
......nm1=n-1 <br> 
......j=1 <br> 
......do 7 i=1,nm1 <br> 
......if(i .ge. j) go to 5 <br> 
......tr=ar(j) <br> 
......ti=ai(j) <br> 
......ar(j)=ar(i) <br> 
......ai(j)=ai(i) <br> 
......ar(i)=tr <br> 
......ai(i)=ti <br> 
....5 k=nv2 <br> 
....6 if(k .ge. j) go to 7 <br> 
......j=j-k <br> 
......k=k/2 <br> 
......go to 6 <br> 
....7 j=j+k <br> 
......pi=3.14159265358979d0 <br> 
......do 20 l=1,mex <br> 
......le=2**l <br> 
......le1=le/2 <br> 
......ur=1.0 <br> 
......ui=0. <br> 
......wr=dcos(pi/le1 ) <br> 
......wi=dsin (pi/le1) <br> 
......do 20 j=1,le1 <br> 
......do 10 i=j,n,le <br> 
......ip=i+le1 <br> 
......tr=ar(ip)*ur - ai(ip)*ui <br> 
......ti=ai(ip)*ur + ar(ip)*ui <br> 
......ar(ip)=ar(i)-tr <br> 
......ai(ip)=ai(i) - ti <br> 
......ar(i)=ar(i)+tr <br> 
......ai(i)=ai(i)+ti <br> 
...10 continue <br> 
......utemp=ur <br> 
......ur=ur*wr - ui*wi <br> 
......ui=ui*wr + utemp*wi <br> 
...20 continue <br> 
......return <br> 
......end <br>