module defpi implicit none real(8),parameter::pi=3.14159265358979323846D0 end module defpi program f90_PHASE !********************************************* ! Fourier Phase Spectrum !********************************************* use defpi implicit none integer::i,k ! integer::ndata !Number of original data integer::nn !Number of calculated data (powers of 2) real(8)::dt !Time increment real(8),allocatable::accinp(:) !Inputted original data real(8),allocatable::xr(:) !Real part of calculated data real(8),allocatable::xi(:) !Imaginary part of calculated data real(8),allocatable::fp(:) !Fourier phase spectrum real(8),allocatable::frq(:) !Fourier phase spectrum real(8)::pwave character::fnameR*50,fnameW*50,fnameP*50 character::strcom*200,dummy*50,fmt1*200 character::linebuff*1000 call getarg(1,fnameR) !Input data file name call getarg(2,fnameW) !Phase angle output file name call getarg(3,fnameP) !Phase wave output file name open(11,file=fnameR,status='old') read(11,'(a)') strcom read(11,*) dummy,dt read(11,*) dummy,ndata !To set number of data ! (nn is equal to 2 to the power of positive integer number ! and nn is greater than or equal to number of data) nn=2 do nn=nn*2 if(ndata*1<=nn)exit end do allocate(accinp(1:nn)) allocate(xr(1:nn)) allocate(xi(1:nn)) allocate(fp(1:nn)) allocate(frq(1:nn/2+1)) do k=1,nn/2+1 !set frequency frq(k)=1.0D0/dt/dble(nn)*dble(k-1) end do do i=1,nn accinp(i)=0.0D0 xr(i)=0.0D0 xi(i)=0.0D0 end do do i=1,ndata read(11,*) accinp(i) xr(i)=accinp(i) end do close(11) call FFT(nn,xr,xi) !Fourier transform do i=1,nn !To divide returned values by number of data xr(i)=xr(i)/dble(nn) xi(i)=xi(i)/dble(nn) fp(i)=0.0D0 if(1.0D-30=k) j=j-k k=k/2 if(k==0)exit end do j=j+k end do end subroutine FFT subroutine del_spaces(s) character (*), intent (inout) :: s character (len=len(s)) tmp integer i, j j = 1 do i = 1, len(s) if (s(i:i)==' ') cycle tmp(j:j) = s(i:i) j = j + 1 end do s = tmp(1:j-1) end subroutine del_spaces end program f90_PHASE