*********************************************************************** * cft512 - 2-dimensional FFT (512x512) * * coordinate of the DC component of the output : (257,257) * * istat : input,integer direction of FFT * * istat = 1 * * a : input,real real part of the original map * * output,real cos of Fourier map * * b : input,real imaginary part of the original map * * output,real sin of Fourier map * * istat = -1 * * a : input,real(512,512) cos of Fourier map * * output,real(512,512) real part of the output map * * b : input,real(512,512) sin of Fourier map * * output,real(512,512) imaginary part of the output map * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine cft512(a,b,istat) c real a(512,512),b(512,512),c(512,512) integer naxes(2) c ***** normal FFT ***** if (istat.eq.1) then c ***** replace image center (257,257) -> (0,0) ***** do 300 il=1,512 do 320 i=1,512 c(i,il)=a(mod(i+255,512)+1,mod(il+255,512)+1) b(i,il)=0. 320 continue 300 continue c ***** FFT ***** naxes(1)=512 naxes(2)=512 m=2 call cft(c,b,naxes,m,1,icon) c write(6,'('' icon ='',i7)') icon c do 350 il=1,512,64 c 350 write(6,'('' '',8f9.2)') (c(i,il)/512./512.,i=1,512,64) c do 360 il=1,512,64 c 360 write(6,'('' '',8f9.2)') (b(i,il)/512./512.,i=1,512,64) c ***** replace image center (0,0) -> (257,257) ***** do 410 il=1,512 do 420 i=1,512 a(i,il)=c(mod(i+255,512)+1,mod(il+255,512)+1)/512./512. 420 continue 410 continue do 430 il=1,512 do 430 i=1,512 430 c(i,il)=b(mod(i+255,512)+1,mod(il+255,512)+1)/512./512. c do 440 il=1,512 do 440 i=1,512 440 b(i,il)=c(i,il) c end if c ***** inverse FFT ***** if (istat.eq.-1) then c ***** replace image center (257,257) -> (0,0) ***** do 100 l=1,512 do 100 i=1,512 100 c(i,l)=a(mod(i+255,512)+1,mod(l+255,512)+1) do 110 l=1,512 do 110 i=1,512 110 a(i,l)=c(i,l) do 120 l=1,512 do 120 i=1,512 120 c(i,l)=b(mod(i+255,512)+1,mod(l+255,512)+1) c ***** inverse FFT ***** m=2 naxes(1)=512 naxes(2)=512 call cft(a,c,naxes,m,-1,icon) c write(6,'('' icon ='',i7)') icon c ***** replace image center (0,0) -> (257,257) ***** do 200 l=1,512 do 200 i=1,512 200 b(i,l)=c(mod(i+255,512)+1,mod(l+255,512)+1) do 210 l=1,512 do 210 i=1,512 210 c(i,l)=a(i,l) do 220 l=1,512 do 220 i=1,512 220 a(i,l)=c(mod(i+255,512)+1,mod(l+255,512)+1) c c do 250 il=1,512,64 c 250 write(6,'('' '',8f9.2)') (a(i,il),i=1,512,64) c do 260 il=1,512,64 c 260 write(6,'('' '',8f9.2)') (b(i,il),i=1,512,64) c end if c return end *********************************************************************** * cft128 - 2-dimensional FFT (128x128) * * coordinate of the DC component of the output : (65,65) * * istat : input,integer direction of FFT * * istat = 1 * * a : input,real real part of the original map * * output,real cos of Fourier map * * b : input,real imaginary part of the original map * * output,real sin of Fourier map * * istat = -1 * * a : input,real(128,128) cos of Fourier map * * output,real(128,128) real part of the output map * * b : input,real(128,128) sin of Fourier map * * output,real(128,128) imaginary part of the output map * * June 22,1992 Y. Hanaoka * * Feb. 26,1999 K. Fujiki * *********************************************************************** c subroutine cft128(a,b,istat) c real a(128,128),b(128,128),c(128,128) integer naxes(2) c ***** normal FFT ***** if (istat.eq.1) then c ***** replace image center (65,65) -> (0,0) ***** do 300 il=1,128 do 320 i=1,128 c(i,il)=a(mod(i+63,128)+1,mod(il+63,128)+1) b(i,il)=0. 320 continue 300 continue c ***** FFT ***** naxes(1)=128 naxes(2)=128 m=2 call cft(c,b,naxes,m,1,icon) c write(6,'('' icon ='',i7)') icon c do 350 il=1,128,64 c 350 write(6,'('' '',8f9.2)') (c(i,il)/128./128.,i=1,128,64) c do 360 il=1,128,64 c 360 write(6,'('' '',8f9.2)') (b(i,il)/128./128.,i=1,128,64) c ***** replace image center (0,0) -> (65,65) ***** do 410 il=1,128 do 420 i=1,128 a(i,il)=c(mod(i+63,128)+1,mod(il+63,128)+1)/128./128. 420 continue 410 continue do 430 il=1,128 do 430 i=1,128 430 c(i,il)=b(mod(i+63,128)+1,mod(il+63,128)+1)/128./128. c do 440 il=1,128 do 440 i=1,128 440 b(i,il)=c(i,il) c end if c ***** inverse FFT ***** if (istat.eq.-1) then c ***** replace image center (65,65) -> (0,0) ***** do 100 l=1,128 do 100 i=1,128 100 c(i,l)=a(mod(i+63,128)+1,mod(l+63,128)+1) do 110 l=1,128 do 110 i=1,128 110 a(i,l)=c(i,l) do 120 l=1,128 do 120 i=1,128 120 c(i,l)=b(mod(i+63,128)+1,mod(l+63,128)+1) c ***** inverse FFT ***** m=2 naxes(1)=128 naxes(2)=128 call cft(a,c,naxes,m,-1,icon) c write(6,'('' icon ='',i7)') icon c ***** replace image center (0,0) -> (65,65) ***** do 200 l=1,128 do 200 i=1,128 200 b(i,l)=c(mod(i+63,128)+1,mod(l+63,128)+1) do 210 l=1,128 do 210 i=1,128 210 c(i,l)=a(i,l) do 220 l=1,128 do 220 i=1,128 220 a(i,l)=c(mod(i+63,128)+1,mod(l+63,128)+1) c c do 250 il=1,128,64 c 250 write(6,'('' '',8f9.2)') (a(i,il),i=1,128,64) c do 260 il=1,128,64 c 260 write(6,'('' '',8f9.2)') (b(i,il),i=1,128,64) c end if c return end *********************************************************************** * cft1024 - 2-dimensional FFT (1024x1024) * * coordinate of the DC component of the output : (513,513) * * istat : input,integer direction of FFT * * istat = 1 * * a : input,real real part of the original map * * output,real cos of Fourier map * * b : input,real imaginary part of the original map * * output,real sin of Fourier map * * istat = -1 * * a : input,real(1024,1024) cos of Fourier map * * output,real(1024,1024) real part of the output map * * b : input,real(1024,1024) sin of Fourier map * * output,real(1024,1024) imaginary part of the output map * * June 22,1992 Y. Hanaoka * * Feb. 26,1999 K. Fujiki * *********************************************************************** c subroutine cft1024(a,b,istat) c real a(1024,1024),b(1024,1024),c(1024,1024) integer naxes(2) c ***** normal FFT ***** if (istat.eq.1) then c ***** replace image center (513,513) -> (0,0) ***** do 300 il=1,1024 do 320 i=1,1024 c(i,il)=a(mod(i+511,1024)+1,mod(il+511,1024)+1) b(i,il)=0. 320 continue 300 continue c ***** FFT ***** naxes(1)=1024 naxes(2)=1024 m=2 call cft(c,b,naxes,m,1,icon) c write(6,'('' icon ='',i7)') icon c do 350 il=1,1024,64 c 350 write(6,'('' '',8f9.2)') (c(i,il)/1024./1024.,i=1,1024,64) c do 360 il=1,1024,64 c 360 write(6,'('' '',8f9.2)') (b(i,il)/1024./1024.,i=1,1024,64) c ***** replace image center (0,0) -> (513,513) ***** do 410 il=1,1024 do 420 i=1,1024 a(i,il)=c(mod(i+511,1024)+1,mod(il+511,1024)+1)/1024./1024. 420 continue 410 continue do 430 il=1,1024 do 430 i=1,1024 430 c(i,il)=b(mod(i+511,1024)+1,mod(il+511,1024)+1)/1024./1024. c do 440 il=1,1024 do 440 i=1,1024 440 b(i,il)=c(i,il) c end if c ***** inverse FFT ***** if (istat.eq.-1) then c ***** replace image center (513,513) -> (0,0) ***** do 100 l=1,1024 do 100 i=1,1024 100 c(i,l)=a(mod(i+511,1024)+1,mod(l+511,1024)+1) do 110 l=1,1024 do 110 i=1,1024 110 a(i,l)=c(i,l) do 120 l=1,1024 do 120 i=1,1024 120 c(i,l)=b(mod(i+511,1024)+1,mod(l+511,1024)+1) c ***** inverse FFT ***** m=2 naxes(1)=1024 naxes(2)=1024 call cft(a,c,naxes,m,-1,icon) c write(6,'('' icon ='',i7)') icon c ***** replace image center (0,0) -> (513,513) ***** do 200 l=1,1024 do 200 i=1,1024 200 b(i,l)=c(mod(i+511,1024)+1,mod(l+511,1024)+1) do 210 l=1,1024 do 210 i=1,1024 210 c(i,l)=a(i,l) do 220 l=1,1024 do 220 i=1,1024 220 a(i,l)=c(mod(i+511,1024)+1,mod(l+511,1024)+1) c c do 250 il=1,1024,64 c 250 write(6,'('' '',8f9.2)') (a(i,il),i=1,1024,64) c do 260 il=1,1024,64 c 260 write(6,'('' '',8f9.2)') (b(i,il),i=1,1024,64) c end if c return end