c======================================================================= c powspec: This routine reads in a times series of acceleration c expectation values and calculates a power spectrum. Two c routines from Numerical Recipes in Fortran are used (realft and c four1). This program is designed to compute spectra for systems c of the form H = H0 + V(t)*x*cos(w0*t+theta). The power spectrum c is produced as a function of w/w0. The maximum value of w/w0 c for which the power spectrum is calculated is given by c c hmax = N/(2*ncyc), where N is the number of points in the c time series and ncyc is the number of cycles of w0 in the c time series. c======================================================================= program powspec implicit none c----------------------------------------------------------------------- c argument parsing variables c----------------------------------------------------------------------- integer iargc, i4arg real*8 r8arg, r8_never parameter (r8_never = -1.0d-60) c----------------------------------------------------------------------- c maxn: maximum number of points in the time series c----------------------------------------------------------------------- integer maxn, n, i, p2 parameter (maxn = 16 384) c----------------------------------------------------------------------- c vectors for storing times, accel, and a modified form of accel c that will be used in the fourier transform c----------------------------------------------------------------------- real*8 data(maxn) real*8 time(maxn), accel(maxn) c----------------------------------------------------------------------- c w, hmax, ncyc as defined above c ho: harmonic order c pow: power at a given harmonic order c----------------------------------------------------------------------- real*8 w, hmax integer ncyc real*8 ho, pow, pi c----------------------------------------------------------------------- c argument parsing: w, ncyc c----------------------------------------------------------------------- if (iargc() .lt. 2) go to 900 w = r8arg(1,r8_never) if (w .eq. r8_never) go to 900 ncyc = i4arg(2, -1) if (ncyc .lt. 0) go to 900 c----------------------------------------------------------------------- c read in the acceleration time series c----------------------------------------------------------------------- call rdacc(time,accel,n,maxn) c----------------------------------------------------------------------- c Error message if there is no data read by rdacc c----------------------------------------------------------------------- if (n .lt. 0) then write(0,*) 'No values found in standard input.' stop end if c----------------------------------------------------------------------- c Check to make sure that n is a power of 2 (maximum of 2**14) c The fourier transform won't work unless it is c----------------------------------------------------------------------- do p2 = 1, 14 if (n .eq. 2**p2) go to 1 end do write(0,*) 'Number of data points must be a power of 2!' write(0,*) 'Unable to compute Fourier Transform.' stop 1 continue c----------------------------------------------------------------------- c Define pi, hmax c----------------------------------------------------------------------- pi = 4.0d0*atan(1.0d0) hmax = dble(n)/(2*dble(ncyc)) c----------------------------------------------------------------------- c Modify accel values for fourier transform c----------------------------------------------------------------------- do i = 1, n data(i) = accel(i)*(1.0d0-((2.0d0*dble(i)-dble(n-1))/ & dble(n+1))**2) end do c----------------------------------------------------------------------- c Calculate fourier transform using Num. Rec. routine 'realft' c----------------------------------------------------------------------- call realft(data,n,1) c----------------------------------------------------------------------- c Compute power and write ho, power to standard out c----------------------------------------------------------------------- pow = pi*data(1)**2/(w**2*hmax**2) write(*,100) 0.0e0, log10(pow) 100 format(1p,2e25.16) do i = 2, n/2 ho = dble(i-1)/dble(ncyc) pow = pi*(data(2*i-1)**2 + data(2*i)**2)/(w**2*hmax**2) write(*,200) ho, log10(pow) 200 format(1p,2e25.16) end do stop 900 continue write(0,*) 'usage: powspec ' write(0,*) write(0,*) 'w: frequency of driving field' write(0,*) 'ncyc: # of cycles in time series' write(0,*) write(0,*) 'reads acceleration data from standard in' write(0,*) 'writes spectrum to standard out' stop end c======================================================================= c History: dvfrom (by Matt Choptuik) c this is a subroutine that reads an acceleration time series from c standard input and puts the times and accleration values into c the vectors time(n) and accel(n). Note that the dimension n c is determined by how many values are in standard input (up to c maxn) c======================================================================= subroutine rdacc(time,accel,n,maxn) implicit none integer n, maxn real*8 time(maxn), accel(maxn) c----------------------------------------------------------------------- c Temporaries for reading values of time and accel c----------------------------------------------------------------------- real*8 tn, an c----------------------------------------------------------------------- c return code from read statement c----------------------------------------------------------------------- integer rc, ustdin parameter (ustdin = 5) c----------------------------------------------------------------------- c initialize c----------------------------------------------------------------------- n = 0 c----------------------------------------------------------------------- c implied do loop to read all values c----------------------------------------------------------------------- 100 continue read(ustdin,*,iostat=rc,end=200) tn, an if (rc .eq. 0) then n = n + 1 c----------------------------------------------------------------------- c Print error message if there are more than maxn points c----------------------------------------------------------------------- if (n .gt. maxn) then write(0,*) 'rdacc: Read maximum of ',maxn, & ' from standard input.' n = maxn return end if c----------------------------------------------------------------------- c Inject read values into vectors c----------------------------------------------------------------------- time(n) = tn accel(n) = an end if go to 100 c----------------------------------------------------------------------- c break out of loop when end-of-file is reached c----------------------------------------------------------------------- 200 continue return end c======================================================================= c This is a double-precision version of the 'realft' routine c in Numerical Recipes in Fortran. It calculates the Fourier c transform of a set of n real-valued data points. The vector c containing the data is replaced by a vector containing the c fourier components. It will also calculate the inverse c transform of a complex data set if it is the transform of c real data. c======================================================================= subroutine realft(data,n,isign) implicit none integer isign, n real*8 data(n) integer i, i1, i2, i3, i4, n2p3 real*8 c1, c2, h1i, h1r, h2i, h2r real*8 theta, wi, wpi, wpr, wr, wtemp isign = 1 c----------------------------------------------------------------------- c Initialize the recurrence c----------------------------------------------------------------------- theta = 3.141592653589793d0/dble(n/2) c1 = 0.5d0 c----------------------------------------------------------------------- c Perform the forward fourier transform c Note: 'four1' performs the transform c----------------------------------------------------------------------- c2 = -0.5d0 call four1(data,n/2,+1) wpr = -2.0d0*sin(0.5d0*theta)**2 wpi = sin(theta) wr = 1.0d0 + wpr wi = wpi n2p3 = n + 3 c----------------------------------------------------------------------- c Case for i = 1 done separately below c----------------------------------------------------------------------- do i = 2, n/4 i1 = 2*i - 1 i2 = i1 + 1 i3 = n2p3 - i2 i4 = i3 + 1 c----------------------------------------------------------------------- c Separates the two transforms from the combined transform c----------------------------------------------------------------------- h1r = c1*(data(i1)+data(i3)) h1i = c1*(data(i2)-data(i4)) h2r = -c2*(data(i2)+data(i4)) h2i = c2*(data(i1)-data(i3)) c----------------------------------------------------------------------- c Separate transforms are recombined to form the true transform c of the original data c----------------------------------------------------------------------- data(i1) = h1r+wr*h2r-wi*h2i data(i2) = h1i+wr*h2i+wi*h2r data(i3) = h1r-wr*h2r+wi*h2i data(i4) = -h1i+wr*h2i+wi*h2r c----------------------------------------------------------------------- c The recurrence c----------------------------------------------------------------------- wtemp = wr wr = wr*wpr-wi*wpi+wr wi = wi*wpr+wtemp*wpi+wi end do c----------------------------------------------------------------------- c Squeeze first and last data c together to get all values within the original array c----------------------------------------------------------------------- h1r = data(1) data(1) = h1r+data(2) data(2) = h1r-data(2) return end c======================================================================= c This routine computes the fourier transform of a complex array c of length nn. It replaces the original data array with the c array of complex fourier coefficients. c c Note: nn must be a power of 2! c c This routine is also from Num. Rec. c======================================================================= subroutine four1(data,nn,isign) implicit none integer isign, nn real*8 data(2*nn) integer i, istep, j, m, mmax, n real*8 tempi, tempr real*8 theta, wi, wpi, wpr, wr, wtemp n = 2*nn j = 1 c----------------------------------------------------------------------- c Perform bit reversal of nn c----------------------------------------------------------------------- do i = 1, n, 2 if (j .gt. i) then tempr = data(j) tempi = data(j+1) data(j) = data(i) data(j+1) = data(i+1) data(i) = tempr data(i+1) = tempi end if m = n/2 1 if ((m .ge. 2) .and. (j .gt. m)) then j = j-m m = m/2 goto 1 end if j = j+m end do c----------------------------------------------------------------------- c Implement the Danielson-Lanczos algorithm. The outer loop is c executed log2(nn) times c----------------------------------------------------------------------- mmax = 2 2 if (n .gt. mmax) then istep = 2*mmax c----------------------------------------------------------------------- c Initialize for trigonometric recurrence c----------------------------------------------------------------------- theta = 6.28318530717959d0/(isign*mmax) wpr = -2.0d0*sin(0.5d0*theta)**2 wpi = sin(theta) wr = 1.0d0 wi = 0.0d0 c----------------------------------------------------------------------- c Two nested inner loops c----------------------------------------------------------------------- do m = 1, mmax, 2 do i = m, n, istep c----------------------------------------------------------------------- c This is the Danielson-Lanczos formula c----------------------------------------------------------------------- j = i + mmax tempr = wr*data(j)-wi*data(j+1) tempi = wr*data(j+1)+wi*data(j) data(j) = data(i)-tempr data(j+1) = data(i+1)-tempi data(i) = data(i)+tempr data(i+1) = data(i+1)+tempi end do c----------------------------------------------------------------------- c trigonometric recurrence c----------------------------------------------------------------------- wtemp = wr wr = wr*wpr-wi*wpi+wr wi = wi*wpr+wtemp*wpi+wi end do mmax = istep goto 2 end if return end