program cross c c Sine multitaper version for cross spectra c c$$$$ calls ascan getdat stats units getint spcmat adapt2 c$$$$ calls getchr clear plot getone yule2 later lag c c c Uses input paradigm for command structure. c One nice feature allows the user to specify the units of c the measurements and then to get a properly normalized c PSD out in correct units. c c Common blocks of ascan: parameter (inmx=200) character *80 input(inmx) common /store/ input common /ndict/ iecho,nin c c Common blocks of getdat, state c Note: to increase the length of series that can be handled c raise mxx in every routine. Memory requirements are: c 11*mxx + constant parameter (mxx=100 000) common /series/ nx,dt,x(mxx,2) common /state/ xbar(2),varx(2),sxmin(2),sxmax(2),ebarx(2) c c Common block from adapt2 parameter (nfmx=mxx/2 + 2) common /result/ nt,nf,fNyq,df,sxy(nfmx,4),kopt(nfmx) c c parameter (deg=180.0/3.1415926535) character*64 unit,spout data unit/' '/, spout/'fort.1'/ data i1,i2,i3/1,1,1/, ktop/15/ data ntimes/3/, initap/2/, lord/0/ c ce dt=1.0 nx=0 c c Read in the commands 1000 call ascan c c c Replot spectra from previous calculations call getone('replot', dum, jagain) if (jagain .eq. 1 .and. nx.gt.0) then call plot goto 1000 endif c Input the data series call getdat call stats c c call units(i1,i2,i3, unit) write(6,'(/(a,2g12.4,1x,a))') $' Mean values of series:',xbar,unit(i1:i2), $' Total variances of the series:',varx,unit(i2:i3), $' RMS signal levels:',sqrt(varx(1)),sqrt(varx(2)), $ unit(i1:i2) write(6,'(a,g12.4,1x,a)') $' Sampling interval of series:',dt, unit(1:i1) c write(6,'(a)')' ','Mean values have been removed from both series' $ c c Prewhiten both series with same AR filter derived from 1st or c 2nd time series according to the sign of lord. Both series x(.,.) c overwritten in yule2. call getint('prewhiten', lord, kord) if (lord .ne. 0) then nser=(3 - sign(1,lord))/2 nord=abs(lord) call yule2(1, nser, nord) write(6,'(/a,i3/a,i2/a)') $ 'Both series prewhitened with an AR filter, order ',nord, $ 'derived from the spectrum of series',nser,' ' endif c c c Decide on the type of calculation initap=3.0 + 0.3*sqrt(real(nx)) call getint('adapt', ntimes, nugat) call getint('tapers', initap, nugat) if (later('tapers', 'adapt') .gt. 0) then c Uniform number of tapers requested write(6,'(a,i6,a)')'All spectra computed with ',initap,' tapers' $ call spcmat(1, -initap) c Adaptive tapering requested - default else mxtap=2*sqrt(real(nx)) c fact=1.0 call getone('smooth', fact, nugat) if (fact .ne. 1.0) write(6,'(a,f5.2)') $ ' Variance weighted by ',fact**2 if (fact.lt.0.1 .or. fact.gt.5.0) write(6,'(a)') $ '>>>> Smoothing factor is probably too large or too small' c write(6,'(a,i6,a)') $ 'Adaptive estimation with ', ntimes,' iterations' initap=initap*sqrt(fact) call adapt2(ntimes, initap, fact) endif c c Restore spectral matrix if either series was 1st differenced or c prewhitened call undiff call yule2(2, nser, nord) c fNyq=0.5/dt df=fNyq/(nf - 1) c c Normalize by variance in x1 = area under psd power=0.5*(sxy(1,1) + sxy(nf,1)) do 1600 k=2, nf-1 power=power + sxy(k,1) 1600 continue const=varx(1)/(power*df) do 1750 k=1, nf do 1700 i=1, 4 sxy(k,i)=sxy(k,i)*const 1700 continue 1750 continue c c c Print results write(6,'(/(a,g12.4,1x,a))') $' Nyquist frequency:',fNyq,unit(1:i1), $' Frequency spacing:',df,unit(1:i1) write(6,'(/(a,i6))') $'Number of terms taken in series:',nt, $' Number of frequencies:',nf c c c c Write results to fort.1 or other designated file call getchr('output', spout, named) if (spout.eq.'-none') spout='off' if (spout .ne. 'off') then if (named .gt. 0) then close (unit=1) open (unit=1, file=spout) call clear('output') endif endif c kmin=nf kmax=0 kbar=0 do 2000 j=1, nf gamsq=(sxy(j,3)**2 + sxy(j,4)**2)/(sxy(j,1)*sxy(j,2)) gain=sqrt(gamsq*sxy(j,2)/sxy(j,1)) phase=deg* atan2( sxy(j,4), sxy(j,3)) if (spout .ne. 'off') $ write(1,'(1p,e11.5,2g11.4,3g12.4,0p,f7.4,f8.2,i5)') $ (j-1)*df,(sxy(j,k),k=1,4),gain,gamsq,phase,kopt(j) kmin=min(kopt(j), kmin) kmax=max(kopt(j), kmax) kbar=kbar + kopt(j) call lag(j, phase, gamsq) 2000 continue kbar=kbar/nf write(6,'(a,i6,a,2(i5,a))') $ ' Average number of tapers:',kbar, $ ' [min',kmin,' max',kmax,']' c call lag(0, dphidf, gamsq) clag = 2.7778e-3/df write(6,'(7x,a,f5.1,a)') $'Series 2 lags series 1 by ',-dphidf*clag,' time steps' c c Plot the results call plot c if (spout.eq. 'off')write(6,'(/a)')'No results written to disk' if (spout.ne. 'off')write(6,'(/(a))') $ 'Results written to '//spout(1:50), $ ' f, S11, S22, Re S12, Im S12, gain, g12**2, PHI, Kopt' goto 1000 end c c__________________________________________________________________ subroutine plot c$$$$ calls units getchr getone getint getarr system clear c c Creates plotxy plotfile of current spectrum, or of all c spectra. Invokes plotxy via 'system' call c c Common blocks of getdat, state c Notice: to increase the length of series that can be treated, c simply raise mxx in every routine. Memory requirements are: c K*mxx + constant parameter (mxx=100 000) common /series/ nx,dt,x(mxx,2) common /state/ xbar(2),varx(2),sxmin(2),sxmax(2),ebarx(2) c c Common block from adapt2 parameter (nfmx=mxx/2 + 2) common /result/ nt,nf,fNyq,df,sxy(nfmx,4),kopt(nfmx) c parameter (deg=180.0/3.1415926535) character*64 unit,splot,title,perfq*12 dimension caylim(2) data unit/' '/, splot/'fort.2'/,title/' '/, nplot/0/ data i1,i2,i3/1,1,1/ data caylim /0.0, 0.0/, lord/0/ c c Get units of the series and frequency ce call units(i1,i2,i3, unit) c c Make a plotfile for results c nplot=nplot+1 call getchr('plot', splot, nugat) if (splot.eq.'-none') splot='off' if (splot .eq. 'off') return c open (unit=2, file=splot) c c c See if a subinterval of freq axis is requested caylim(1)=-100.0 call getarr('detail', caylim, 2, ndet) call getarr('replot', caylim, 2, ndet) call clear('replot') nf1= caylim(1)/df + 1.0 nf2= caylim(2)/df - 1.0 if (.not.(1.le.nf1.and.nf1.lt.nf2.and.nf2.le.nf))then nf1=1 nf2=nf endif c Ascertain if log frequency is desired in plotted output call getint('logfreq', nugat, logf) if (logf .ge. 0) nf1=max(nf1, 2) c c c Write plotfile to unit 2 c c Power spectra to lower panel of graph write(2,'(a)')'title \\com\\','logxy 2','frame','file *', $ 'char 0.08','ylim 2.5','xlim 4' if (logf .ge. 0) write(2,'(a)') 'logxy' if (unit(1:1) .eq. ' ' .and. dt.eq.1.0) then write(2,'(a)') 'xlabel Frequency in 1/\\D\\t' elseif (unit(1:1) .eq. ' ') then write(2, '(a)') 'xlabel Frequency ' elseif (unit(1:2) .eq. 'Hz') then write(2, '(a)') 'xlabel Frequency Hz' else write(2, '(3a)') 'xlabel Frequency ',unit(2:i1-1),'\\sup{-1}' endif if (i2 .le. i1+1) then write(2,'(a)') 'ylabel PSD in arbitrary units' else perfq=unit(2:i1) if (unit(1:i1) .eq. 'Hz') perfq='/Hz' write(2,'(4a)') $ 'ylabel PSD ',unit(i1:i2-1),'\\sup{2}',perfq endif c Note if series are prewhitened call getint('prewhiten', lord, kord) if (lord .ne. 0) write(2,'(a)') $'note (0.2 3.6i) \\0008\\'// $'Both series prewhitened' c Write out the spectra write(2,'(a/a,3g14.6)') 'mode 1', $ 'affine ',df,(nf1-2)*df,1.0,0.0 write(2,'(a,i6/(1p,7g11.4))') 'read',nf2-nf1+1, $ (sxy(k,1), k=nf1, nf2) write(2,'(a/a,i6/(1p,7g11.4))') 'color dark red','read', $ nf2-nf1+1, (sxy(k,2), k=nf1, nf2) write(2,'(a)') 'color black','plot 1.5 5.0','note' c c Plot the gain write(2,'(a/a,i6)') 'color blue','read',nf2-nf1+1 do 1110 j=nf1,nf2 gain=sqrt(sxy(j,3)**2 + sxy(j,4)**2)/sxy(j,1) write(2,'(e11.4)') gain 1110 continue write(2,'(a)')'ylim 2.5 0 0','ylab Gain G\\sub{12}(f)', $'frame -xaxis right off', $'xlab','plot 0 0','color black','frame on +xaxis left' c c c Upper panel of graph: phase and coherence/gain call getchr('title', title, ititi) if (ititi.lt. 0) call getchr('file', title, ignor) write(2,'(a,a)')'title ',title c c First plot confidence value in gray; fudge of 1.5 for parabolic c weighting in specmat p=0.95 call getone('confidence', p, nugat) if (0.0 .lt. p .and. p.lt. 1.0) $ write(2,'(a/a/a,i6/(10f8.4))') $ 'color dark gray','logxy linlin','read,',nf2-nf1+1, $ (1.0-(1.0-p)**(1.5/(kopt(j)-1)), j=nf1,nf2) if (logf .ge. 0) write(2,'(a)') 'logxy loglin' c c Now the coherence write(2,'(a/a,i6)') 'color black','read',nf2-nf1+1 do 1100 j=nf1,nf2 gamsq=(sxy(j,3)**2 + sxy(j,4)**2)/(sxy(j,1)*sxy(j,2)) write(2,'(f9.5)') gamsq 1100 continue write(2,'(a)')'ylim 2.5 0 1','ylab \\g\\\\sup{2}', $'xlab','note','plot' c c Now the phase call getone('unwrap', dum, nrap) write(2,'(a/a,i6)')'color dark green','read',nf2-nf1+1 do 1200 j=nf1,nf2 phase=deg* atan2( sxy(j,4), sxy(j,3)) if (nrap .ge. 0 .and. j.gt.1) then dph=phase-oldph if (dph .gt. 185.0) phase=phase-360.0 if (dph .lt.-185.0) phase=phase+360.0 endif oldph=phase if (j.eq.1 .and. nrap.lt.0) phase=360.0 write(2,'(f9.2)') phase 1200 continue write(2,'(a)') 'frame off top right -xaxis','ylim 2.5 -185 185 45' if (nrap .ge. 0) write(2,'(a)')'ylim 2.5 0 0 45' c c Finish off write(2,'(a)') 'ylab Phase','plot 0 0','stop' close (unit=2) c write(6,'(/a)')'Plotxy file written to '//splot c c c WARNING: 'system' is a SUNOS special command not generally available c in Fortran which starts a Unix shell. Here it runs plotxy. If no c equivalent is available, delete and run the plot command from c another window or after cross has terminated. call getone('hold', dum, ihold) if (ihold .ge. 0) return call system('plotxy < '//splot//' > /dev/null ; '// $' gv mypost ') c c return end c______________________________________________________________________ subroutine undiff c$$$$ calls getint c Corrects spectral matrix if one of the series has been 1st c differenced before analysis as a form of primitive prewhitening. c parameter (mxx=100 000) parameter (nfmx=mxx/2 + 2) parameter (pi=3.141592653589793d0) complex s12,ghat common /result/ nt,nf,fNyq,df,sxy(nfmx,4),kopt(nfmx) c c ipik=0 call getint('difference', ipik, npik) if (npik .lt. 0) return ipik=max(1, min(2, ipik)) c del=sign(pi/(2.0*nf), 1.5-ipik) do 1500 j=2, nf ghat=2.0*(1.0 - exp(cmplx(0.0, j*del))) s12=cmplx(sxy(j,3), sxy(j,4))/ghat sxy(j,3)=real(s12) sxy(j,4)=aimag(s12) sxy(j,ipik)=sxy(j,ipik)/abs(ghat)**2 1500 continue c Copy 2nd freq est into 1st do 1600 k=1, 4 sxy(1,k)=sxy(2,k) 1600 continue return end c_______________________________________________________________________ subroutine lag(j, phase, gamsq) c$$$$ calls nothing c If j=0 return mean and std deviation in phase, gamsq c if j=1 initialize c if j>1 find various sums related to the sum of delta-phi c weighted by coherence c common /hold/ pho,gmo,sw,swx c c Estimates mean and its variance if (j .eq. 0) then phase = swx/sw c Initializes sums and other parameters elseif (j .eq. 1) then pho = phase gmo = gamsq sw = 0.0 swx = 0.0 else x = phase - pho c Exclude phase wrap if (abs(x) .le. 200.0) then w = gamsq + gmo sw = sw + w swx = swx + w*x endif pho = phase gmo = gamsq endif return end c_______________________________________________________________________ subroutine yule2(iphase, nser, nord) c$$$$ calls toepl slowft nearn fft c Prewhitening of the two spectra. See Percival and Walden, c "Spectral Analaysis ...", Chap 9, 1993. c c nord = number of AR coeffs to be found. c nser=1, find pilot spectral estimates on series 1 for c autocovariances ar() for Yule-Walker equations. c nser = 2, use the second series. c iphase=1: Calculate AR filter and prewhiten data series c iphase=2; correct spectra for AR filtering c c Double-length FT of must be ready in y() in /pool/ c Both series x(,1), x(,2) are prewhitened with the same AR filter. parameter (pi=3.141592653589793) parameter (mxx=100 000) parameter (nfmx=mxx/2 + 2) complex y,zz common /series/ nx,dt,x(mxx,2) common /result/ nt,nf,fNyq,df,sxy(nfmx,4),kopt(nfmx) common /pool/ y(2*mxx+2,2) parameter (ndmx=20) dimension row(2*ndmx),rhs(ndmx),ar(ndmx),wrk(2*ndmx) dimension yre(4*mxx+4,2) equivalence (y(1,1), yre(1,1)) save ar c ce if (nord .le. 1) return c if (iphase .eq. 2) goto 2000 c iphase=1: Prewhiten the 2 time series nrd=min(ndmx-1, nord) if (nrd .ne. nord) write(6,'(a,i3)') $'>>> Order of prewhitening filter reduced to ',nrd c c c Make pilot spectral estimate for series nser c c Perform double-length complex FFT on data and save in y(.,1). ny=2*nearn(nx/2) n2=2*ny nf=ny/2 + 1 do 1100 j=1, ny y(j,1)=x(j,nser) 1100 continue do 1150 j=ny+1, n2 y(j,1)=0.0 1150 continue c c Perform mixed-radix FFTs call fft(yre(1,1), yre(2,1), n2,n2,n2, 2, ierr) if (ierr .ne. 0) then write(6,'(a,i7)')'>>> FFT routine fails for n =',n2, $ ' This error should never arise' endif c c c Make sine multitaper spectral estimate, stored in sxy(.,1) klim=nint(5.0 + 0.2*sqrt(real(nx))) c Loop over frequency do 1250 m=0, nf-1 j=m+1 m2=2*m sxy(j,1)=0.0 c c Loop over the number of tapers do 1200 k=1, klim j1=mod(m2+n2-k, n2) j2=mod(m2+k, n2) zz=y(j1+1,1) - y(j2+1,1) sxy(j,1)=sxy(j,1) + real(zz)**2 + aimag(zz)**2 1200 continue sxy(j,1)=sxy(j,1)/real(klim) 1250 continue c c Prewhitening nff=nf-1 c Fourier transform PSD estimate to find autocovariance function. c Load values into Toeplitz matrix and rhs vector row(nrd+2)=0.0 do 1300 i=nrd+1, 1 , -1 omi=pi*(i-1)/nff call slowft(nff, sxy(1,1), omi, row(i), st) row(nrd+i)=row(i) rhs(i)=row(i+1) 1300 continue c c Solve symmetric Toeplitz system for filter coeffs ar() ar(1)=-1.0 call toepl(nrd, row, rhs, wrk, ar(2), iok) c c Convolve filter ar() into both series x(), overwriting do 1600 k=1, 2 do 1500 j=nx, 1, -1 sum=x(j,k) do 1400 jf=2, nrd+1 jx=max(1, j-jf+1) sum=sum - ar(jf)*x(jx,k) 1400 continue x(j,k)=sum 1500 continue 1600 continue c return c c iphase=2: Restore spectra 2000 dom=pi/nf c Divide out from sx() the AR spectrum of filter coeffs ar() do 2100 k=1, nf call slowft(nrd+1, ar, (k-1)*dom, ct,st) fark=1.0 / (ct**2 + st**2) sxy(k,1)=sxy(k,1) * fark sxy(k,2)=sxy(k,2) * fark sxy(k,3)=sxy(k,3) * fark sxy(k,4)=sxy(k,4) * fark 2100 continue return end c________________________________________________________________ subroutine toepl(n, r, y, ff, x, iok) c$$$$ calls nothing c Solves linear system in Toeplitz form. c n order of matrix. c r(1) ... r(n) 1st row of Toeplitz array. Overwritten. c r(n+1) ... r(2*n) 1st col of Toeplitz array. Overwritten. c y() the rhs vector. c ff() working array of size at least 2*n. c x() the solution vector. c iok 1 if everything ok, 0 if singular principal minor. c dimension r(*),x(*),y(*),ff(2,*) c c Re-order r into crazy pattern: c r(1),r(2), ... r(n): the first row of the matrix BACKWARDS! c r(n),r(n+1), ... r(2*n-1): the first column of the matrix. ce do 1100 i=1, n/2 ri=r(i) r(i)=r(n-i+1) r(n-i+1)=ri 1100 continue do 1150 i=n+1, 2*n r(i-1)=r(i) 1150 continue c iok=1 if(r(n).eq.0.) goto 99 x(1)=y(1)/r(n) if(n.eq.1)return ff(1,1)=r(n-1)/r(n) ff(2,1)=r(n+1)/r(n) do 15 m=1,n m1=m+1 sxn=-y(m1) sd=-r(n) do 11 j=1,m sxn=sxn+r(n+m1-j)*x(j) sd=sd+r(n+m1-j)*ff(1,m-j+1) 11 continue if(sd.eq.0.)goto 99 x(m1)=sxn/sd do 12 j=1,m x(j)=x(j)-x(m1)*ff(1,m-j+1) 12 continue if(m1.eq.n) return sgn=-r(n-m1) shn=-r(n+m1) sgd=-r(n) do 13 j=1,m sgn=sgn+r(n+j-m1)*ff(1,j) shn=shn+r(n+m1-j)*ff(2,j) sgd=sgd+r(n+j-m1)*ff(2,m-j+1) 13 continue if(sd.eq.0..or.sgd.eq.0.)goto 99 ff(1,m1)=sgn/sgd ff(2,m1)=shn/sd k=m m2=(m+1)/2 pp=ff(1,m1) qq=ff(2,m1) do 14 j=1,m2 pt1=ff(1,j) pt2=ff(1,k) qt1=ff(2,j) qt2=ff(2,k) ff(1,j)=pt1-pp*qt2 ff(1,k)=pt2-pp*qt1 ff(2,j)=qt1-qq*pt2 ff(2,k)=qt2-qq*pt1 k=k-1 14 continue 15 continue c 99 iok=0 return end c_______________________________________________________________________ subroutine slowft(n, x, om, ct, st) c$$$$ calls nothing c c Calculates Fourier transform of real sequence x(i),i=1,...n c at angular frequency om normalized so that nyquist=pi: c cmplx(ct, st) = sum {k=0 to n-1} x(k+1)*exp(i*k*om) c The sine transform is returned in st and the cosine transform in c ct. Algorithm is that of Goertzal with modifications by Gentleman, c Comp.J. 1969. The transform is not normalized. To normalize c one-sided ft, divide by sqrt(data length). For positive om, the ft c is defined as ct-(0.,1.)st or like slatec cfftf c c Slightly modified from Slatec library c implicit double precision (a-h,o-z) real x, om, ct, st parameter (pi=3.141592653589793238d0,tp=2.d0*pi) dimension x(n) ce omega=om np1=n+1 l=6.d0*omega/tp s=sin(omega) a=0.d0 c=0.d0 d=0.d0 e=0.d0 if(l.eq.0)then c Recursion for low frequencies (.lt. nyq/3) b=-4.d0*sin(omega/2.d0)**2 do 10 k=1,n c=a d=e a=x(np1-k)+b*d+c 10 e=a+d elseif(l.eq.1)then c Regular Goertzal algorithm for intermediate frequencies b=2.d0*cos(omega) do 20 k=1,n a=x(np1-k)+b*e-d d=e 20 e=a else c Recursion for high frequencies (.gt. 2*nyq/3) b=4.d0*cos(omega/2.d0)**2 do 30 k=1,n c=a d=e a=x(np1-k)+b*d-c 30 e=a-d endif st=-s*d ct=a-b*d/2.d0 return end c_______________________________________________________________________ c======================================================================= c Unit D: Get the Data c======================================================================= subroutine getdat c$$$$ calls getarr getchr getone getint spline eval ljust akima evlak c c Gets the file name and reads data into x(.,2) c May read from 2 columns in a single file or from two files c Spline interpolates if time unevenly spaced. c character*64 name, temp, even, oldev, spln*3 parameter (mxx=100 000) dimension tab(40),column(3),skp(2) common /series/ nx,dt,x(mxx,2) common /state/ xbar(2),varx(2),sxmin(2),sxmax(2),ebarx(2) common /pool/ y(8*mxx+8) c Memory for splines shared with other arrays dimension t(mxx),u(mxx,2),ddu(mxx),wrk(mxx) equivalence (t(1),y(1)),(ddu(1),y(mxx+1)), $(u(1,1),y(2*mxx+1)),(wrk(1),y(4*mxx+1)) c data name(1:1)/' '/,nterm/mxx/, column/1,2,3/ data spln/'nat'/, oldev/' '/ save nterm ce call getchr('file', name, nln) if (nln .eq. 0) then write(6,'(a)')'>>> Program cannot continue without a file name' stop endif kbl=index(name(1:nln), ' ') nfiles=min(2, kbl+1) if (kbl .eq. 0) kbl=nln c c If there are two files, squeeze out excess blanks in name if (nfiles .eq. 2) then call ljust(name(kbl:nln), temp) name(kbl+1:nln)=temp write(6,'(a)') $ 'Data will be read from two separate files:',name column(2)=1 endif i1=1 i2=kbl c c Determine number of records to be skipped before reading skp(1)=0.0 call getarr('skip', skp, nfiles, none) if (none .le. 1) skp(2)=skp(1) c c do 1800 nfl=1, nfiles c open(unit=11, file=name(i1:i2), status='OLD', err=3100) i1=kbl+1 i2=nln c iskip=skp(nfl) do 1010 j=1, iskip read (11,'(a1)',end=3600) dum 1010 continue if (none.eq.2 .or. iskip.gt.0) write(6,'(a,i7)') $ 'Number of records skipped before reading data:',iskip c c Is series to be decimated? call getint('decimate', idec, idc) if (idc .le. 0) idec=1 idec=max(1, idec) if (idec .gt. 1 .and. nfiles .gt.1) goto 3700 if (idec .gt. 1) then wsum=0.0 do 1020 j=1, idec wrk(j)=sin(3.1415926*j/(idec+1.0))**1.5 wsum=wsum + wrk(j) 1020 continue do 1022 j=1, idec wrk(j)=wrk(j)/wsum 1022 continue endif c c Get number of terms to be read nt=mxx*idec call getint('nterms', nterm, nte) if (nte .gt. 0) nt=min(mxx*idec, nterm) c c Get the sampling interval call getone('interval', dt, intdt) dt = dt*idec if (idec .gt. 1) write(6,'(a,g15.4)') $ 'After decimation new sampling interval:',dt c c Get columns to be read call getarr('column', column, 3, kols) c c Check for spline interpolation and conditions call getchr('spline', spln, ispln) c if (ispln .ge. 0 .and. intdt .le. 0) goto 3200 if (ispln .ge. 0 .and. kols.gt.0 .and. kols.lt.3) goto 3210 if (ispln .ge. 0 .and. nfiles .eq. 2) goto 3220 if (ispln .ge. 0 .and. idec .gt. 1) goto 3710 if (kols .eq. 1) goto 3300 c kol1=nint(column(1)) kol2=nint(column(2)) koln=max(kol1, kol2) c Check if columns OK to read if (.not. (koln .ge. 1 .and. koln .le. 40 .and. $ min(kol1,kol2).ge.1)) goto 3500 c c c c Both series read from the same file if (nfiles .eq. 1 .and. ispln.lt.0 .and. idec.le.1) then write(6,'(a,2i4)')'Series read from columns ',kol1,kol2 do 1050 j=1, nt read (11,*,err=3000,end=1700) (tab(l),l=1,koln) x(j,1)=tab(kol1) x(j,2)=tab(kol2) 1050 continue c c Both series from 1 file and decimation requested elseif (nfiles .eq. 1 .and. ispln.lt.0 .and. idec.gt.1) then write(6,'(a,2i4/a,i4)')'Series read from columns ',kol1,kol2, $ 'Both series decimated by',idec, $ 'All parameters listed for series AFTER decimation' np=nt/idec do 1090 j=1, np x(j,1)=0.0 x(j,2)=0.0 do 1060 k=1, idec read (11,*,err=3000,end=1095) (tab(l),l=1,koln) x(j,1)=wrk(k)*tab(kol1) + x(j,1) x(j,2)=wrk(k)*tab(kol2) + x(j,2) 1060 continue 1090 continue j=np+1 1095 nx=j-1 goto 1900 c c c Both series from same file and spline interpolation requested elseif (ispln .ge. 0) then kol3=nint(column(3)) koln=max(koln, kol3) write(6,'(a,3i4)') $ 'Time and data series read from columns',kol1,kol2,kol3 do 1100 j=1, nt read (11,*,err=3000, end=1110) (tab(l),l=1,koln) t(j)=tab(kol1) u(j,1)=tab(kol2) u(j,2)=tab(kol3) if (j.ge.2 .and. t(j).le.t(j-1)) goto 3400 1100 continue 1110 nj=j-1 close (unit=11) write(6,'(a,i7)')'Number of terms before interpolation:',nj do 1300 ks=1, 2 c Select natural or Akima splines, and interpolate if (spln .eq. 'nat') then if (ks.eq.1)write(6,*)'Natural spline interpolation' call spline(nj, t, u(1,ks), ddu, wrk) do 1200 j=1, mxx tx=t(1) + (j-1)*dt if (tx .gt. t(nj)) goto 1210 x(j,ks)=eval(tx, nj, t, u(1,ks), ddu) 1200 continue else if (ks.eq.1)write(6,*)'Akima splines used to interpolate' call akima(nj, t, u(1,ks), wrk) do 1205 j=1, mxx tx=t(1) + (j-1)*dt if (tx .gt. t(nj)) goto 1210 x(j,ks)=evlak(tx, nj, t, u(1,ks), wrk) 1205 continue endif j=mxx+1 1210 nx=j-1 1300 continue goto 1800 c c c Series read from two files else koln=nint(column(nfl)) write(6,'(2(a,i3))')'Reading from data file number',nfl, $ ', column ',koln do 1500 j=1, nt read (11,*, err=3000, end=1700) (tab(l),l=1,koln) x(j,nfl)=tab(koln) 1500 continue endif c j=nt+1 1700 nx=j-1 if (nfl .eq. 1) nfirst=nx if (nx .ne. nfirst .and. nfl.eq. 2) write(6,'(a,2i6)') $ '>>> Different number of terms in the two files:',nfirst,nx, $ ' The longer series is truncated to parity' nx=min(nx, nfirst) close(unit=11) 1800 continue c c If requested write interpolated series to a diskfile call getchr('validate', even, nval) if (nval .gt. 0) then if (oldev .ne. even) then if (oldev(1:1) .ne. ' ') close (unit=4) open (unit=4, file=even) oldev=even endif write(4,'(3g13.5)') (t(1)+(j-1)*dt,x(j,1),x(j,2),j=1,nx) write(6,'(/a,a)') $ 'Interpolated series written to: ',even(1:nval) elseif (nval .eq. 0) then write(6,'(a)')'>>> A filename is needed with validate' endif c c c 1900 write(6,'(/a,i6)')'Number of terms in data series:',nx,' ' if (nx .eq. mxx) write(6,'(a)') $'>>> Array space filled: series probably truncated' c c 1st difference both series if requested - primitive prewhitening ipik=0 call getint('difference', ipik, idiff) if (idiff .lt. 0) return ipik=max(1, min(2, ipik)) do 2200 j=1, nx-1 x(j,ipik)=x(j+1,ipik) - x(j,ipik) 2200 continue nx=nx - 1 write(6,'(a,i3,a)') $'Series',ipik,' is prewhitened by 1st differencing' return c c Error out 3000 write(6,'(a,i7)')'>>> Unreadable number in data file at point',j stop 3100 write(6,'(2a)')'>>> Unable to open the file ',name(i1:i2) stop 3200 write(6,'(a)')'>>> You must set the interval to use spline' stop 3210 write(6,'(a)')'>>> You must provide 3 columns to use spline' stop 3220 write(6,'(a)')'>>> For spline both series must be in 1 file' stop 3300 write(6,'(a)')'>>> You must provide 2 numbers for column' stop 3400 write(6,'(a,i7)')'>>> Time variable not increasing at term:',j stop 3500 write(6,'(a,2g12.4)')'>>> Column numbers unacceptable:', $ column(1),column(2) stop 3600 write(6,'(a)') $'>>> End-of-file encountered during initial skipping' stop 3700 write(6,'(a)')'>>> For decimation both series must be in 1 file' stop 3710 write(6,'(a)')'>>> You may not decimate and spline a series' stop end c_______________________________________________________________________ subroutine stats c$$$$ calls getone c Finds some statistics of the series - removes mean and trend c Saves answers in /state/ parameter (mxx=100 000) common /series/ nx,dt,x(mxx,2) common /state/ xbar(2),varx(2),sxmin(2),sxmax(2),ebarx(2) double precision sum,xmean data detr/-1/ ce c Step through each series in turn do 1900 k=1, 2 sum=0.0 c Find mean do 1100 j=1, nx sum=x(j,k) + sum 1100 continue xmean=sum/nx xbar(k)=xmean c Find variance sum=0.0 do 1500 j=1, nx sum=sum + (x(j,k)-xmean)**2 1500 continue varx(k)=sum/nx c c Removes mean from both series do 1600 j=1, nx x(j,k)=x(j,k) - xmean 1600 continue c c Estimates slope of LS line tx=0.0 tbar=0.5+ 0.5*nx do 1700 j=1, nx tx=tx + (j-tbar)*x(j,k) 1700 continue xn=nx slope=12.0*tx/(xn*(xn**2+6.0*xn-4.0)) write(6,'(a,g12.4,a,i2)') 'Total offset produced by trend: ', $ slope*nx,' in series ',k c Remove the slope if requested call getone('detrend', detr, mark) if (mark .ge. 0)then sum=0.0 do 1800 j=1, nx x(j,k)=x(j,k) - (j-tbar)*slope sum=sum + x(j,k)**2 1800 continue varx(k)=sum/nx if (k.eq.2)write(6,'(a)')'Both series have been detrended' endif 1900 continue * WRITE(3,'(2g12.4)') (x(j,1),x(j,2),j=1,nx) c return end c_______________________________________________________________________ subroutine units(i1,i2,i3, unit) c$$$$ calls getchr c Reads in the names of the units of time variable and signal variable. c On return unit(1:i1) holds time variable unit + a space c unit(i1+1,i2) holds signal variable unit name + space c unit(i2+1,i3) holds signal-squared unit name + space character*64 unit ce call getchr('units', unit(2:64), inut) if (inut .eq. 0) then unit=' ' elseif (inut .ge. 1) then unit(1:1)='/' i1=index(unit, ' ') i2=index(unit(i1+1:64), ' ')+i1 i3=2*i2 - i1 + 2 if (i2.gt.i1+1) unit(i2+1:i3)=unit(i1+1:i2-1)//'^2' if (unit(1:i1) .eq. '/s ') unit(1:i1)='Hz ' endif return end c_______________________________________________________________________ subroutine akima(n, x,f, df) c$$$$ calls nothing c Given the array of samples of a function f and sample points x, c assumed to be montone, generates slopes for an interpolating rule c according to Akima's algorithm (Lancaster and Salkauskas, c Curve and surface fitting, 1986 academic press, p 82). c dimension x(n),f(n),df(n),S(4) c q(u1,x1,u2,x2)=(u1/x1**2-u2/x2**2)/(1.0/x1-1.0/x2) c ce S(1)=(f(2) - f(1))/(x(2) - x(1)) S(2)=S(1) S(3)=S(1) S(4)=(f(3) - f(2))/(x(3) - x(2)) Sn=(f(n) - f(n-1))/(x(n) - x(n-1)) eps=1.0e-6*abs(x(n) - x(1)) do 1200 i=1, n D1=abs(S(2) - S(1)) D2=abs(S(4) - S(3)) df(i)=(D2*S(2) + D1*S(3))/(D1 + D2 + eps) c S(1)=S(2) S(2)=S(3) S(3)=S(4) if (i+3 .le. n) then S(4)=(f(i+3) - f(i+2))/(x(i+3)-x(i+2)) else S(4)=Sn endif 1200 continue c if (n .eq. 2) return c If 3 or more points use gradient from a parabola for 1st & last df df(1)=q(f(2)-f(1),x(2)-x(1),f(3)-f(1),x(3)-x(1)) df(n)=q(f(n-1)-f(n),x(n-1)-x(n),f(n-2)-f(n),x(n-2)-x(n)) return end c____________________________________________________________________ function evlak(y, n,x,f,df) c$$$$ calls nothing c Given a montone increasing array of x values (Note: routine does not c check this), samples of values of the function f and 1st derivative c df at the xs, interpolates by cubic polynomial. The array df can be c found by calling subroutine akima. c c If y falls outside the interval, the polynomial on the c nearest interval is used to extrapolate. dimension x(n), f(n), df(n) c c Search for proper interval initiated at previous call if possible. save init data init/1/ c c Locate sample interval containing y: after y lies in [x(init), c x(init+1)), unless y lies outside [x(1), x(n)] when the interval c is the intervals containing the apprpriate end point. ce init=min(init, n) if (y .gt. x(init)) then do 1100 k=init, n if (x(k) .gt. y) then init=k-1 goto 1300 endif 1100 continue init=n-1 else do 1200 k=init, 1, -1 if (x(k) .le. y) then init=k goto 1300 endif 1200 continue init=1 endif 1300 dx=x(init+1) - x(init) c c Evaluate the cubic interpolator t=(y - x(init))/dx s=1.0 - t evlak=s**2*((1.0 + 2.0*t)*f(init) + t*dx*df(init)) + $ t**2*((1.0 + 2.0*s)*f(init+1) - s*dx*df(init+1)) return end c____________________________________________________________________ subroutine spcmat(new, ktop) c$$$$ calls fft nearn c Computes the Spectral matrix Sij from the two time series c x(*,2) in /series/, using the same number of sine tapers c at each frequency, kopt(). Based on PSD estimation of c Riedel & Sidorenko, IEEE Tr Sig Pr,43, 188, 1995 c new =0 means no need to fft the series - these are already saved c in yre(,). Otherwise, perform the ffts c ktop = maximum of tapers permitted in estimates c parameter (mxx=100 000) parameter (pi=3.141592653589793d0) parameter (nfmx=mxx/2 + 2) complex y,z1,z2 common /series/ nx,dt,x(mxx,2) common /result/ nt,nf,fNyq,df,sxy(nfmx,4),kopt(nfmx) common /pool/ y(2*mxx+2,2) dimension yre(4*mxx+4,2) equivalence (y(1,1), yre(1,1)) c c c ktop if > 0, largest number in array kopt c if < 0, -ktop = constant value to be used; insert c -ktop into the array kopt. c kopt() nf-vector giving number of tapers to be averaged c into each estimate a frequency number j, if ktop > 0. c Works on: c x(*,2) two colum n-vectors of the time series in /series/ c Returns into /results/: c nf number of frequencies in spectrum c sxy(*,4) Spectral matrix of cross spectra with frequency: c sxy(*,1), sxy(*,2) PSDs of 2 series; c sxy(*,3), sxy(*,4) co and quad spectra of the series. c kopt() is filled if ktop < 0. c c c Up to double the length of the time series ce ny=2*nearn(nx/2) n2=2*ny nf=ny/2 + 1 nt=ny fNyq=0.5/dt df=fNyq/(nf - 1) c c Perform double-length complex FFT on data and save in y. if (new .ge. 1) then do 1100 j=1, ny y(j,1)=x(j,1) y(j,2)=x(j,2) 1100 continue do 1150 j=ny+1, n2 y(j,1)=0.0 y(j,2)=0.0 1150 continue c c Perform mixed-radix FFTs call fft(yre(1,1), yre(2,1), n2,n2,n2, 2, ierr) call fft(yre(1,2), yre(2,2), n2,n2,n2, 2, ierr) if (ierr .ne. 0) then write(6,'(a,i7)')'>>> FFT routine fails for n =',n2, $ ' This error should never arise' endif endif c c Fill kopt vector if necessary if (ktop .lt. 0) then do 1200 j=1, nf kopt(j)=-ktop 1200 continue endif c c Loop over frequency do 1500 m=0, nf-1 j=m+1 klim=kopt(j) ck=1.0/klim**2 m2=2*m c do 1300 i=1, 4 sxy(j,i)=0.0 1300 continue c c Average over the tapers with parabolic weight wk do 1400 k=1, klim j1=mod(m2+n2-k, n2) j2=mod(m2+k, n2) z1=y(j1+1,1) - y(j2+1,1) z2=y(j1+1,2) - y(j2+1,2) wk=1.0 - ck*(k-1)**2 c sxy(j,1)=sxy(j,1) + $ wk*(real(z1)**2 + aimag(z1)**2) sxy(j,2)=sxy(j,2) + $ wk*(real(z2)**2 + aimag(z2)**2) sxy(j,3)=sxy(j,3) + $ wk*(real(z1)*real(z2) + aimag(z1)*aimag(z2)) sxy(j,4)=sxy(j,4) + $ wk*(real(z2)*aimag(z1) - real(z1)*aimag(z2)) 1400 continue * wt=1.0/real(klim) c Exact normalization for parabolic weighting wt=6.0*klim/real(4*klim**2+3*klim-1) do 1450 i=1, 4 sxy(j,i)=sxy(j,i)*wt 1450 continue c 1500 continue c return end c_______________________________________________________________________ subroutine adapt2(ntimes, ntap, fact) c$$$$ calls spcmat north calmer c c Perform adaptive spectral estimation by determing the best c number of tapers, kopt, at each frequency. c c ntimes number of iterations c ntap initial number of tapers c fact log S is smoothed (1=normal, >1 smoother) via chisq misfit c Returns in /result/: c nf number of frequencies in PSD c sx() nf-vector of PSD estimates c kopt() integer nf-vector of number of tapers at each freq c parameter (mxx=100 000) parameter (nfmx=mxx/2 + 2) common /series/ nx,dt,x(mxx,2) common /result/ nt,nf,fNyq,df,sxy(nfmx,4),kopt(nfmx) parameter (my=mxx/2+2) dimension opt(nfmx,2),y(my) c c1, c2=(20*sqrt(1.2))**0.4 are constants for parabolic weighting c in subroutine spcmat; for uniform weighting c1=1, c2=12.0**0.4=2.702 data c1/1.2000/, c2/3.437/ c c Get pilot estimate of spectrum ce call spcmat(1, -ntap) c c do 1600 iter=1, ntimes do 1500 ipsd=1, 2 do 1100 j=1, nf y(j)= log(sxy(j, ipsd)) 1100 continue c do 1300 j=1, nf ispan = kopt(j)*1.4 call north(nf, j-ispan, j+ispan, y, d1y, ddy) R = (ddy + d1y**2)/df**2 ak=kopt(j)/real(2*ispan) phi=720.0*ak**5*(1.0 - 1.286*ak + 0.476*ak**3 - $ 0.0909*ak**5) sigR= sqrt(phi/real(kopt(j))**5) / df**2 optj =c2/(df**4 *( R**2 + 1.4*sigR**2) /fact**2)** 0.2 * write (20,'(1p,7g11.4)') real(ipsd), * $ (j-1)*df,opt(j,ipsd),optj opt(j,ipsd)=optj c 1300 continue c 1500 continue c c Condition the number of tapers call calmer(nf, opt, kopt) c * write(21,'(1p,4g11.4)') * $ ((j-1)*df,opt(j,1),opt(j,2),real(kopt(j)), j=1,nf) c c Recompute spectrum with optimal variable taper numbers call spcmat(0, 1) 1600 continue return end c_______________________________________________________________ subroutine north(n, i1, i2, s, ds, dds) c$$$$ calls nothing c Performs LS fit to s by c a degree-two polynomial in an orthogonal basis. c Returns ds = estimate of 1st derivative ds/dn at center of record c Returns dds = estimate of 2nd derivative c dimension s(*) c L = i2 - i1 + 1 el=L gamma = (el**2 - 1.0)/12.0 u0sq = el u1sq = el*(el**2 - 1.0)/12.0 u2sq = (el*(el**2 - 1.0)*(el**2- 4.0))/180.0 amid= 0.5*(el + 1.0) dot0=0.0 dot1=0.0 dot2=0.0 do 1100 kk=1, L i=kk + i1 - 1 c Negative or excessive index uses even function assumption if (i.le. 0) i=2 - i if (i.gt. n) i=2*n - i dot0 = dot0 + s(i) dot1 = dot1 + (kk - amid) * s(i) dot2 = dot2 + ((kk - amid)**2 - gamma)*s(i) 1100 continue ds = dot1/u1sq dds = 2.0*dot2/u2sq c return end c_______________________________________________________________________ subroutine calmer(nf, opt, kopt) c$$$$ calls curb c Suppresses strong peaks in Kopt curve that are (almost) c always due to S" going through zero. Restricts rate of c growth of weight series to |dK/dn| < 1 parameter (mxx=100 000) parameter (nfmx=mxx/2 + 2) dimension opt(nfmx,2),kopt(nf) c do 1100 j=1,nf opt(j,1)=min(opt(j,1), opt(j,2)) 1100 continue c call curb(nf,opt(1,1)) do 1200 j=1, nf kopt(j)=max(opt(j,1),3.0) 1200 continue return end c_______________________________________________________________________ subroutine curb(n, v) c$$$$ calls nothing c Rewrites the input n-vector v() so that all points lie below c the piece-wise linear function v(k) + abs(j-k), where v(k) c is a local minimum in the original v. c Effectively clips strong peaks and keeps slopes under 1 in c magnitude. dimension v(n) c do 1500 j=2, n-1 c Scan series for local minimum if (v(j).lt.v(j+1) .and. v(j).lt.v(j-1)) then vloc=v(j) c Revise series accordingly do 1200 k=1, n v(k)=min(v(k), vloc+abs(j-k)) 1200 continue endif 1500 continue return end c_______________________________________________________________ function nearn(n) c$$$$ calls nothing c Finds the largest integer less than or equal to n whose c prime factors are in the set {2, 3, 5}. c parameter (c=0.693147181,b=1.098612289,a=1.609437912) ce p=log(real(n)+0.5) ka=p/a kb=p/b kc=p/c c dmin=p do 1500 ja=0, ka do 1400 jb=0, kb do 1300 jc=0, kc d=p -ja*a -jb*b -jc*c if (d .lt. 0.0) goto 1400 if (d .lt. dmin) dmin=d 1300 continue 1400 continue 1500 continue c nearn=nint(exp(p-dmin)) return end c_______________________________________________________________ c======================================================================= c UNIT K: DECODING ROUTINES FOR COMMAND KIT c======================================================================= subroutine ascan c$$$$ calls clear getchr icheck c Input routine for commands c c Reads from the standard input until eof or code "execute " or "quit". c Saves lines in the input store /store/ for later retrieval by c getarr, getone, getint and getchr. c c Prints a glossary upon request c parameter (inmx=200) character *80 input(inmx),line, code*4 common /store/ input common /ndict/ iecho,nin c if (nin .eq. 0) write(6,'(a)') ' ', $'Enter commands for spectral analysis (? for help)' c do 1500 l=nin+1, inmx read(5,'(80a)', end=3000) line if (line(1:4).eq.'cont' .or. line(1:4).eq.'exec') goto 2000 c List a glossary of codes if (line(1:1) .eq. '?') then write(6, '(a/a/(2x,a))') $'Enter commands from the following list:', $'?: Remind me of the command list again', $'execute: Quit reading commands and begin current calculations', $'file filename: Enter filename of time series (mandatory)', $'output filename: Provide the disk filename for results;', $' filename = off means no output', $'plot filename: Provide filename for plotxy file', $'interval dt: Sampling interval of data', $'units time signal: Measurement units of independent variable', $' and of the signal (eg: sec volts)', $'nterms n: Number of data points to be read from file', $'skip k: Skip over k lines before reading data', $'column k1 k2 [k3]: Read series 1 from column k1 and series 2', $' from colum k2; k3 used with spline', $'spline: kind Data unevenly spaced in time; interpolate to', $' interval dt from file with 3 cols: t, x1,x2; kind=nat or Akima' write(6, '(2x,a)') $'validate filename: Write interpolated series to a file', $'detrend: If present, series are detrended', $'review: Display current command stack', $'detail f1 f2: Plot spectra for f1 < freq < f2', $'replot f1 f2: Replot spectra from the previous calculation', $'tapers nt: Number of tapers used in uniform spectral estimate', $'adapt ntimes: Apply adaptive process ntimes', $'smooth s: Increase s>1, or decrease s<1, spectral smoothness', $'unwrap: Remove large phase jumps on the plots', $'confidence p: Plot in gray p-level confidence coherence', $' is not zero', $'prewhiten nw: Prewhiten both series with nw-term AR filter', $' based on 1st series (nw>0) or 2nd (nw<0)', $'title text: Add a title to the plotted information', $'clear command: Delete earlier occurrences of named command', $' ' c c Review the command stack ce elseif (line(1:4) .eq. 'revi') then write(6,'(5x,a)')' ', '=================== ', $ (input(i)(1:60),i=1,nin),'=================== ' c elseif (line(1:1) .ne. ' ') then if (line(1:4) .eq. 'echo') iecho=-iecho c Homonyms translated if (line(1:4) .eq. 'cols') line(1:4)='colu' if (line(1:4) .eq. 'data') line(1:4)='file' if (line(1:4) .eq. 'quit') then write(6,'(/a)') 'Normal termination' stop endif nin=nin + 1 call icheck(line) input(nin)=line c Clear a command and clear clear itself if (line(1:4) .eq. 'clea') then call getchr('clear', code, nugat) call clear (code) nin=nin - 1 endif endif 1500 continue c 2000 continue n1=max(1, nin-24) write(6,'(5x,a)')' ', '=================== ', $(input(i)(1:60),i=n1,nin),'=================== ' return c 3000 stop end c______________________________________________________________________ subroutine icheck(line) c$$$$ calls nothing c Checks the command line against a list; appends a c warning if com is not present in the list. character*80 line, com*4 ce com=line(1:4) if (0 .eq. $+index('adap clea colu conf deta detr diff file inte',com) $+index('maxt nter outp plot prew repl skip spli tape titl',com) $+index('unit unwr vali smoot hold logf deci',com) $) line=line(1:20)//' <<<<<<< Unrecognized command' return end c______________________________________________________________________ subroutine getarr(code, values, nwant, nfound) c$$$$ calls ljust c c Extracts an array of numbers from the input store. That store is c a large array in common /store/ which has been filled earlier. c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned for numbers. c values the real output array of values found. c nwant the maximum number of numbers expected . c nfound the number of numbers actually found in the input. c If the line contains fewer than nwant values, this is the c value returned in nfound. If an error is discovered c nfound=-n, where n is the number of numbers properly c decoded. If there are no numbers after the codeword c nfound=0. Finally, if the code is absent from the store c nfound=-99 and the array is left undisturbed. c parameter (inmx=200) character *80 input(inmx),line,local,char, code*4 common /store/ input common /ndict/ iecho,nin dimension values(*) c c Read the store in reverse order (Thus last entry is obeyed) ce do 1010 lin=nin, 1, -1 line=input(lin) c Check for code if (code .eq. line(1:4)) goto 1020 1010 continue c Code word not found nfound=-99 return c 1020 if (iecho .ge. 1) write(6,'(2a)')'==> ',line n1=index(line, ' ')+1 n2=index(line, '%') n2=80 + min(n2,1)*(n2 - 81) char=line(n1:n2) c do 1500 n=1, nwant call ljust(char, local) lbl=index(local, ' ') if (lbl .eq. 1) goto 1510 read (local, *, err=2000) values(n) char=local(lbl:80) 1500 continue n=nwant+1 1510 nfound=n - 1 return c 2000 write(6,'(a)')' ', $'>>> Unreadable numbers in this input line:',line nfound=1 - n return end c______________________________________________________________ subroutine getone(code, value, nfound) c$$$$ calls getarr c c Extracts a single number from the input store. That store is c a large array in common /store/ which has been filled earlier. c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned for numbers. c value the real output variable containing the desired number. c nfound is 1 if a number is successfully read in; it is zero c the number is absent or unreadable. nfound = -99 if the c code is absent from the input store. c character*4 code dimension v(1) ce call getarr(code, v, 1, nfound) if (nfound .eq. 1) value=v(1) return end c______________________________________________________________ subroutine getint(code, number, nfound) c$$$$ calls getarr c c Extracts a single integer from the input store. c See getone for details. c character*4 code dimension v(1) ce call getarr(code, v, 1, nfound) if (nfound .eq. 1) then number=nint(v(1)) if (number .ne. v(1)) write(6,'(4a,1p,g16.8/a,i10)') $ '>>> Warning: ',code,' expects an integer argument, but', $ ' found: ',v(1),' Returns: ',number,' ' endif return end c______________________________________________________________ subroutine getchr(code, char, nbytes) c$$$$ calls ljust c Extracts a single character variable from the input store. That c store is a large array in common /store/ which has been filled c earlier. c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned. c char the character output variable containing the desired string. c nbytes is the length of the string read; it is zero if c the line is blank after the code. nbytes = -99 if the c code is absent from the input store. c parameter (inmx=200) character *80 input(inmx),line, char*(*), code*4 common /store/ input common /ndict/ iecho,nin c Inspect the store in reverse order (thus reading latest entry) ce do 1010 lin=nin, 1, -1 line=input(lin) c Check for code word if (code .eq. line(1:4)) goto 1020 1010 continue c Code word not found nbytes=-99 return c 1020 if (iecho .ge. 1) write(6,'(2a)')'==> ',line c n1=index(line, ' ')+1 n2=index(line, '%') n2=80 + min(n2,1)*(n2 - 81) c Check for blank string nbytes=0 if (line(n1:n2) .eq. ' ') return c Save in char everything from 1st non-blank to last non-blank call ljust(line(n1:n2), char) nn=min(len(char), n2-n1+1) do 1200 k=nn, 1, -1 if (char(k:k) .ne. ' ') goto 1210 1200 continue 1210 nbytes= k return end c______________________________________________________________ subroutine clear(code) c$$$$ calls nothing c c Removes the last command entry identified by code. c code A 4-byte identifying code. Only lines in the input store c beginning with this code are scanned. c parameter (inmx=200) character *80 input(inmx),line, code*4 common /store/ input common /ndict/ iecho,nin c c Inspect the store in normal order ce do 1010 lin= nin, 1, -1 line=input(lin) c Check code and clear the line if present throughout stack if (code .eq. line(1:4)) then input(lin)=' --' * return endif 1010 continue return end c______________________________________________________________ subroutine ljust(str1, str2) c c$$$$ calls nothing c Left justify: input character string str1 (up to 80 bytes in length) c is left justified, removing leading blanks, and returned in str2. c str1 and str2 can be set to the same variable in the call. character*(*) str1, str2 character*80 local c str2=str1 if (str1 .eq. ' ') return c ls=min(len(str1), len(str2), 80) str2=' ' local(1:ls)=str1 do 1100 j=1, ls if (local(j:j) .ne. ' ') goto 1200 1100 continue 1200 str2(1:ls)=local(j:ls) return end c______________________________________________________________ blockdata iounit common /ndict/ iecho,nin c Don't echo command lines when read data iecho/-1/ c Initial command line count data nin/0/ end c______________________________________________________________________ function later(code1, code2) c$$$$ calls nothing c Returns the difference in order number in the stack of the c most recent occurrence of the commands code1, code2. If the c command isn't present assigns it the order number zero. c Thus later > 0 if code1 occurs after code 2, < 0 if before, c and later=0 if both are absent and code1 and code2 are different. parameter (inmx=200) character *80 input(inmx),line, code1*4,code2*4 common /store/ input common /ndict/ iecho,nin c Inspect the store ce kode1=0 kode2=0 do 1010 lin=1, nin line=input(lin) if (code1 .eq. line(1:4)) kode1=lin if (code2 .eq. line(1:4)) kode2=lin 1010 continue c later=kode1 - kode2 return end c______________________________________________________________ c======================================================================= c Unit F: Fast Fourier transforms of real series, length = 2**k c======================================================================= subroutine realtr(a, b, n, isn) c$$$$$ calls nothing c Rearranges array for fast Fourier transformation of a real series c in place. c Examples c c Example (1) : forward transform of a real series, length n. c call ffttwo(x, n/2) c call realtr(x, x(2), n/2, 2) c Here ffttwo expects complex array with real,imag parts alternating. c result appears as cos,sin,cos,sin... in x(1),x(2),...x(n+1),x(n+2) c note that x must be dimensioned at least n+2 to hold the trans c -forms of the Nyquist frequency. c c Example (2) : inverse transform back to a real series. c Complex Fourier transform in x is stored as it would be created c by using realtr to perform forward transform. c call realtr(x, x(2), n/2, -2) c call fftinv(x, n/2) c The routine is not restricted to values of n that are powers of 2. c Also abs(isn) is the increment between real and imag parts in the c array x. Thus if x is declared complex in most compilers isn=+-2. c If x is to be stored as two separate real arrays abs(isn) is the c distance between real(x(1)) and aimag(x(1)). c dimension a(*),b(*) parameter (piby2=3.1415926535/2) real im ce inc = abs(isn) nk = n * inc + 2 nh = nk / 2 sd = piby2 / n cd = 2. * sin(sd)**2 sd = sin(sd+sd) sn = 0. if (isn .gt. 0) then cn = 1. a(nk-1) = a(1) b(nk-1) = b(1) else cn = -1. sd = -sd endif do 30 j=1,nh,inc k = nk - j aa = a(j) + a(k) ab = a(j) - a(k) ba = b(j) + b(k) bb = b(j) - b(k) re = cn * ba + sn * ab im = sn * ba - cn * ab b(k) = im - bb b(j) = im + bb a(k) = aa - re a(j) = aa + re aa = cn - (cd * cn + sd * sn) sn = (sd * cn - cd * sn) + sn cn = aa 30 continue c return end c_______________________________________________________________________ subroutine ffttwo(a, n) c$$$$ calls nothing c Routine for finding discrete Fourier transform of a complex series. c Modified from Cooley et al IEEE Trans of Education 12 no 1 (1969) c Argument list c a is a complex array containing the series to be transformed. The c results are returned in a over-writing the original series. c n is the number of terms in the series. It must be a power of 2 or c the routine will print a message and return without doing anything. c Replaces a(m+1), for m=0,1, ... n-1 with c sum from k=0 to n-1 [a(k+1)*exp(-2*pi*(0,1)*k*m/n)] c c Note indices are displaced with respect of frequency, e.g. zero c frequency is held in a(1). c dimension a(n) complex a,u,w,t data pi/3.1415926535/ ce m=log(real(n))/0.693147 if (n .ne. 2**m) go to 1000 nv2=n/2 nm1=n-1 j=1 do 8 i=1, nm1 if (i.lt.j) then t=a(j) a(j)=a(i) a(i)=t endif k=nv2 6 if (k.ge.j) go to 7 j=j-k k=k/2 go to 6 7 j=j+k 8 continue c le=1 do 30 l=1, m le1=le le=2*le u=(1., 0.) b=pi/le1 w=cmplx(cos(b), sin(b)) do 20 j=1, le1 do 10 i=j, n, le ip=i+le1 t=a(ip)*u a(ip)=a(i)-t a(i) =a(i)+t 10 continue u=u*w 20 continue 30 continue return 1000 write(6,'(a,i7,a/(a))')'>>> ', n,' is not a power of 2:', $'FFTTWO cannot transform a series of this length.', $'Series returned unchanged' return end c_______________________________________________________________________ subroutine fft(a, b, ntot, n, nspan, isn, ierr) c Multivariate complex Fourier transform, computed in place using c mixed-radix fast Fourier transform algorithm. By R. C. Singleton, c Stanford Research Institute, Sept. 1968. Arrays a and b originally c hold the real and imaginary components of the data, and return the c real and imaginary components of the resulting Fourier c coefficients. Multivariate data is indexed according to the Fortran c array element successor function, without limit on the number of c implied multiple subscripts. The subroutine is called once for each c variate. The calls for a multivariate transform may be in any c order. c ntot is the total number of complex data values. c n is the dimension of the current variable. c nspan/n is the spacing of consecutive data values while indexing c the current variable. c The integer ierr is an error return indicator. It is normally zero, c but is set to 1 if the number of terms cannot be factored in the c space available. If it is permissible the appropriate action at this c stage is to enter fft again after having reduced the length of the c series by one term. The sign of isn determines the sign of the c complex exponential, and the magnitude of isn is normally one. A c tri-variate transform with a(n1,n2,n3), b(n1,n2,n3) is computed by c call fft(a, b, n1*n2*n3, n1, n1, 1) c call fft(a, b, n1*n2*n3, n2, n1*n2, 1) c call fft(a, b, n1*n2*n3, n3, n1*n2*n3, 1) c For a single-variate c transform, ntot = n = nspan = (number of complex data values), c e.g. call fft(a, b, n, n, n, 1) c With most Fortran compilers the data c can alternatively be stored in a single complex array a, then the c magnitude of isn changed to 2 to give the correct indexing c increment and a(2) used to pass the initial address for the sequence c of imaginary values, e.g. c call fft(a, a(2), ntot, n, nspan, 2) c Arrays at(maxf), ck(maxf), c bt(maxf), sk(maxf), and np(maxp) are used for temporary storage. c if the available storage is insufficient, the program is terminated c by the error return option c maxf must be .ge. the maximum prime factor of n. maxp must be .gt. c the number of prime factors of n. In addition, if the square-free c portion k of n has two or more prime factors, then maxp must be .ge. c k-1. c Array storage in nfac for a maximum of 15 prime factors of n. c If n has more than one square-free factor, the product of the c square-free factors must be .le. 210 parameter (mmaxf=23, maxp=209) c Array storage for maximum prime factor of 23 dimension at(mmaxf),ck(mmaxf),bt(mmaxf),sk(mmaxf) dimension nfac(11),np(maxp) dimension a(*),b(*) equivalence (i,ii) c maxf=mmaxf ierr=0 if(n .lt. 2) return inc=isn c72=0.30901699437494742 s72=0.95105651629515357 s120=0.86602540378443865 rad=6.2831853071796 if(isn .ge. 0) goto 10 s72=-s72 s120=-s120 rad=-rad inc=-inc 10 nt=inc*ntot ks=inc*nspan kspan=ks nn=nt-inc jc=ks/n radf=rad*float(jc)*0.5 i=0 jf=0 c determine the factors of n m=0 k=n goto 20 15 m=m+1 nfac(m)=4 k=k/16 20 if(k-(k/16)*16 .eq. 0) goto 15 j=3 jj=9 goto 30 25 m=m+1 nfac(m)=j k=k/jj 30 if(mod(k,jj) .eq. 0) goto 25 j=j+2 jj=j**2 if(jj .le. k) goto 30 if(k .gt. 4) goto 40 kt=m nfac(m+1)=k if(k .ne. 1) m=m+1 goto 80 40 if(k-(k/4)*4 .ne. 0) goto 50 m=m+1 nfac(m)=2 k=k/4 50 kt=m j=2 60 if(mod(k,j) .ne. 0) goto 70 m=m+1 nfac(m)=j k=k/j 70 j=((j+1)/2)*2+1 if(j .le. k) goto 60 80 if(kt .eq. 0) goto 100 j=kt 90 m=m+1 nfac(m)=nfac(j) j=j-1 if(j .ne. 0) goto 90 c compute Fourier transform 100 sd=radf/float(kspan) cd=2.0*sin(sd)**2 sd=sin(sd+sd) kk=1 i=i+1 if(nfac(i) .ne. 2) goto 400 c transform for factor of 2 (including rotation factor) kspan=kspan/2 k1=kspan+2 210 k2=kk+kspan ak=a(k2) bk=b(k2) a(k2)=a(kk)-ak b(k2)=b(kk)-bk a(kk)=a(kk)+ak b(kk)=b(kk)+bk kk=k2+kspan if(kk .le. nn) goto 210 kk=kk-nn if(kk .le. jc) goto 210 if(kk .gt. kspan) goto 800 220 c1=1.0-cd s1=sd 230 k2=kk+kspan ak=a(kk)-a(k2) bk=b(kk)-b(k2) a(kk)=a(kk)+a(k2) b(kk)=b(kk)+b(k2) a(k2)=c1*ak-s1*bk b(k2)=s1*ak+c1*bk kk=k2+kspan if(kk .lt. nt) goto 230 k2=kk-nt c1=-c1 kk=k1-k2 if(kk .gt. k2) goto 230 ak=cd*c1+sd*s1 s1=(sd*c1-cd*s1)+s1 c1=c1-ak kk=kk+jc if(kk .lt. k2) goto 230 k1=k1+inc+inc kk=(k1-kspan)/2+jc if(kk .le. jc+jc) goto 220 goto 100 c transform for factor of 3 (optional code) 320 k1=kk+kspan k2=k1+kspan ak=a(kk) bk=b(kk) aj=a(k1)+a(k2) bj=b(k1)+b(k2) a(kk)=ak+aj b(kk)=bk+bj ak=-0.5*aj+ak bk=-0.5*bj+bk aj=(a(k1)-a(k2))*s120 bj=(b(k1)-b(k2))*s120 a(k1)=ak-bj b(k1)=bk+aj a(k2)=ak+bj b(k2)=bk-aj kk=k2+kspan if(kk .lt. nn) goto 320 kk=kk-nn if(kk .le. kspan) goto 320 goto 700 c transform for factor of 4 400 if(nfac(i) .ne. 4) goto 600 kspnn=kspan kspan=kspan/4 410 c1=1.0 s1=0 420 k1=kk+kspan k2=k1+kspan k3=k2+kspan akp=a(kk)+a(k2) akm=a(kk)-a(k2) ajp=a(k1)+a(k3) ajm=a(k1)-a(k3) a(kk)=akp+ajp ajp=akp-ajp bkp=b(kk)+b(k2) bkm=b(kk)-b(k2) bjp=b(k1)+b(k3) bjm=b(k1)-b(k3) b(kk)=bkp+bjp bjp=bkp-bjp if(isn .lt. 0) goto 450 akp=akm-bjm akm=akm+bjm bkp=bkm+ajm bkm=bkm-ajm if(s1 .eq. 0) goto 460 430 a(k1)=akp*c1-bkp*s1 b(k1)=akp*s1+bkp*c1 a(k2)=ajp*c2-bjp*s2 b(k2)=ajp*s2+bjp*c2 a(k3)=akm*c3-bkm*s3 b(k3)=akm*s3+bkm*c3 kk=k3+kspan if(kk .le. nt) goto 420 440 c2=cd*c1+sd*s1 s1=(sd*c1-cd*s1)+s1 c1=c1-c2 c2=c1**2-s1**2 s2=2.0*c1*s1 c3=c2*c1-s2*s1 s3=c2*s1+s2*c1 kk=kk-nt+jc if(kk .le. kspan) goto 420 kk=kk-kspan+inc if(kk .le. jc) goto 410 if(kspan .eq. jc) goto 800 goto 100 450 akp=akm+bjm akm=akm-bjm bkp=bkm-ajm bkm=bkm+ajm if(s1 .ne. 0) goto 430 460 a(k1)=akp b(k1)=bkp a(k2)=ajp b(k2)=bjp a(k3)=akm b(k3)=bkm kk=k3+kspan if(kk .le. nt) goto 420 goto 440 c transform for factor of 5 (optional code) 510 c2=c72**2-s72**2 s2=2.0*c72*s72 520 k1=kk+kspan k2=k1+kspan k3=k2+kspan k4=k3+kspan akp=a(k1)+a(k4) akm=a(k1)-a(k4) bkp=b(k1)+b(k4) bkm=b(k1)-b(k4) ajp=a(k2)+a(k3) ajm=a(k2)-a(k3) bjp=b(k2)+b(k3) bjm=b(k2)-b(k3) aa=a(kk) bb=b(kk) a(kk)=aa+akp+ajp b(kk)=bb+bkp+bjp ak=akp*c72+ajp*c2+aa bk=bkp*c72+bjp*c2+bb aj=akm*s72+ajm*s2 bj=bkm*s72+bjm*s2 a(k1)=ak-bj a(k4)=ak+bj b(k1)=bk+aj b(k4)=bk-aj ak=akp*c2+ajp*c72+aa bk=bkp*c2+bjp*c72+bb aj=akm*s2-ajm*s72 bj=bkm*s2-bjm*s72 a(k2)=ak-bj a(k3)=ak+bj b(k2)=bk+aj b(k3)=bk-aj kk=k4+kspan if(kk .lt. nn) goto 520 kk=kk-nn if(kk .le. kspan) goto 520 goto 700 c transform for odd factors 600 k=nfac(i) kspnn=kspan kspan=kspan/k if(k .eq. 3) goto 320 if(k .eq. 5) goto 510 if(k .eq. jf) goto 640 jf=k s1=rad/float(k) c1=cos(s1) s1=sin(s1) if(jf .gt. maxf) goto 998 ck(jf)=1.0 sk(jf)=0.0 j=1 630 ck(j)=ck(k)*c1+sk(k)*s1 sk(j)=ck(k)*s1-sk(k)*c1 k=k-1 ck(k)=ck(j) sk(k)=-sk(j) j=j+1 if(j .lt. k) goto 630 640 k1=kk k2=kk+kspnn aa=a(kk) bb=b(kk) ak=aa bk=bb j=1 k1=k1+kspan 650 k2=k2-kspan j=j+1 at(j)=a(k1)+a(k2) ak=at(j)+ak bt(j)=b(k1)+b(k2) bk=bt(j)+bk j=j+1 at(j)=a(k1)-a(k2) bt(j)=b(k1)-b(k2) k1=k1+kspan if(k1 .lt. k2) goto 650 a(kk)=ak b(kk)=bk k1=kk k2=kk+kspnn j=1 660 k1=k1+kspan k2=k2-kspan jj=j ak=aa bk=bb aj=0.0 bj=0.0 k=1 670 k=k+1 ak=at(k)*ck(jj)+ak bk=bt(k)*ck(jj)+bk k=k+1 aj=at(k)*sk(jj)+aj bj=bt(k)*sk(jj)+bj jj=jj+j if(jj .gt. jf) jj=jj-jf if(k .lt. jf) goto 670 k=jf-j a(k1)=ak-bj b(k1)=bk+aj a(k2)=ak+bj b(k2)=bk-aj j=j+1 if(j .lt. k) goto 660 kk=kk+kspnn if(kk .le. nn) goto 640 kk=kk-nn if(kk .le. kspan) goto 640 c multiply by rotation factor (except for factors of 2 and 4) 700 if(i .eq. m) goto 800 kk=jc+1 710 c2=1.0-cd s1=sd 720 c1=c2 s2=s1 kk=kk+kspan 730 ak=a(kk) a(kk)=c2*ak-s2*b(kk) b(kk)=s2*ak+c2*b(kk) kk=kk+kspnn if(kk .le. nt) goto 730 ak=s1*s2 s2=s1*c2+c1*s2 c2=c1*c2-ak kk=kk-nt+kspan if(kk .le. kspnn) goto 730 c2=c1-(cd*c1+sd*s1) s1=s1+(sd*c1-cd*s1) kk=kk-kspnn+jc if(kk .le. kspan) goto 720 kk=kk-kspan+jc+inc if(kk .le. jc+jc) goto 710 goto 100 c permute the results to normal order---done in two stages c permutation for square factors of n 800 np(1)=ks if(kt .eq. 0) goto 890 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) goto 810 k3=np(k+1) kspan=np(2) kk=jc+1 k2=kspan+1 j=1 if(n .ne. ntot) goto 850 c permutation for single-variate transform (optional code) 820 ak=a(kk) a(kk)=a(k2) a(k2)=ak bk=b(kk) b(kk)=b(k2) b(k2)=bk kk=kk+inc k2=kspan+k2 if(k2 .lt. ks) goto 820 830 k2=k2-np(j) j=j+1 k2=np(j+1)+k2 if(k2 .gt. np(j)) goto 830 j=1 840 if(kk .lt. k2) goto 820 kk=kk+inc k2=kspan+k2 if(k2 .lt. ks) goto 840 if(kk .lt. ks) goto 830 jc=k3 goto 890 c permutation for multivariate transform 850 k=kk+jc 860 ak=a(kk) a(kk)=a(k2) a(k2)=ak bk=b(kk) b(kk)=b(k2) b(k2)=bk kk=kk+inc k2=k2+inc if(kk .lt. k) goto 860 kk=kk+ks-jc k2=k2+ks-jc if(kk .lt. nt) goto 850 k2=k2-nt+kspan kk=kk-nt+jc if(k2 .lt. ks) goto 850 870 k2=k2-np(j) j=j+1 k2=np(j+1)+k2 if(k2 .gt. np(j)) goto 870 j=1 880 if(kk .lt. k2) goto 850 kk=kk+jc k2=kspan+k2 if(k2 .lt. ks) goto 880 if(kk .lt. ks) goto 870 jc=k3 890 if(2*kt+1 .ge. m) return kspnn=np(kt+1) c permutation for square-free factors of n j=m-kt nfac(j+1)=1 900 nfac(j)=nfac(j)*nfac(j+1) j=j-1 if(j .ne. kt) goto 900 kt=kt+1 nn=nfac(kt)-1 if(nn .gt. maxp) goto 998 jj=0 j=0 goto 906 902 jj=jj-k2 k2=kk k=k+1 kk=nfac(k) 904 jj=kk+jj if(jj .ge. k2) goto 902 np(j)=jj 906 k2=nfac(kt) k=kt+1 kk=nfac(k) j=j+1 if(j .le. nn) goto 904 c determine the permutation cycles of length greater than 1 j=0 goto 914 910 k=kk kk=np(k) np(k)=-kk if(kk .ne. j) goto 910 k3=kk 914 j=j+1 kk=np(j) if(kk .lt. 0) goto 914 if(kk .ne. j) goto 910 np(j)=-j if(j .ne. nn) goto 914 maxf=inc*maxf c reorder a and b, following the permutation cycles goto 950 924 j=j-1 if(np(j) .lt. 0) goto 924 jj=jc 926 kspan=jj if(jj .gt. maxf) kspan=maxf jj=jj-kspan k=np(j) kk=jc*k+ii+jj k1=kk+kspan k2=0 928 k2=k2+1 at(k2)=a(k1) bt(k2)=b(k1) k1=k1-inc if(k1 .ne. kk) goto 928 932 k1=kk+kspan k2=k1-jc*(k+np(k)) k=-np(k) 936 a(k1)=a(k2) b(k1)=b(k2) k1=k1-inc k2=k2-inc if(k1 .ne. kk) goto 936 kk=k2 if(k .ne. j) goto 932 k1=kk+kspan k2=0 940 k2=k2+1 a(k1)=at(k2) b(k1)=bt(k2) k1=k1-inc if(k1 .ne. kk) goto 940 if(jj .ne. 0) goto 926 if(j .ne. 1) goto 924 950 j=k3+1 nt=nt-kspnn ii=nt-inc+1 if(nt .ge. 0) goto 924 return c Error finish, insufficient array storage 998 write(6,*)' Array bounds exceeded within subroutine fft' ierr=1 end c_______________________________________________________________________ c======================================================================= c Unit S: Splines and smoothing splines c======================================================================= subroutine spline(nn, x, u, s, a) c$$$$ calls nothing c Finds array s for spline interpolator eval. c Based on an alogorithm by Richard F. Thompson in 'Spline c interpolation on a digital computer' (x-692-70-261) 1970 Goddard c Space Flight Center. Note these splines are not 'natural' since the c 2nd derivative is not zero x(1) and x(n). c nn number of data points supplied (may be negative, see below) c x() nn-array of x-coordinates where function is sampled. The c sequence xx(1),xx(2)... must be a strictly increasing sequence c but THIS IS NOT CHECKED by spline. c u() nn-vector of function values that are to be interpolated. c s() nn-vector output of 2nd derivative at sample points. c a() nn-vector, working space. c c If the user wishes derivatives at the ends of series to assume c specified values, du(1)/dx, du(nn)/dx must be placed in s1, s2 c and in the call nn=-number of terms in series. Normally a parabola c is fitted through the 1st and last 3 points to find the slopes. c c If 2 points are supplied, a straight lines is fitted. c If 3 points are supplied, a parabola is fitted. c c At evaluation time, straight-line extrapolation is provided c outside the interval (x(1), x(n)). common /startx/ istart dimension x(*),u(*),s(*),a(*) ce q(u1, x1, u2, x2)=(u1/x1**2-u2/x2**2)/(1.0/x1-1.0/x2) c c Assign derivatives at the ends. istart=0 n=iabs(nn) c Series too short for cubic spline - use straight lines. if (n.le. 2) then s(1)=0.0 s(2)=0.0 return endif c q1=q(u(2)-u(1), x(2)-x(1), u(3)-u(1), x(3)-x(1)) qn=q(u(n-1)-u(n), x(n-1)-x(n), u(n-2)-u(n), x(n-2)-x(n)) c Three points - fit a parabola if (n .eq. 3) then s(1)=(qn - q1)/(x(3) - x(1)) s(2)=s(1) s(3)=s(1) return endif c Force the derivatives to take user-defined values if (nn.le.0) then q1=s(1) qn=s(2) endif s(1)=6.0*((u(2)-u(1))/(x(2)-x(1)) - q1) n1= n - 1 do 2000 i=2,n1 s(i)= (u(i-1)/(x(i)-x(i-1)) - u(i)*(1.0/(x(i)-x(i-1))+ + 1.0/(x(i+1)-x(i))) + u(i+1)/(x(i+1)-x(i)))*6.0 2000 continue s(n)=6.0*(qn + (u(n1)-u(n))/(x(n)-x(n1))) a(1)=2.0*(x(2)-x(1)) a(2)=1.5*(x(2)-x(1)) + 2.0*(x(3)-x(2)) s(2)=s(2) - 0.5*s(1) do 3000 i=3,n1 c=(x(i)-x(i-1))/a(i-1) a(i)=2.0*(x(i+1)-x(i-1)) - c*(x(i)-x(i-1)) s(i)=s(i) - c*s(i-1) 3000 continue c=(x(n)-x(n1))/a(n1) a(n)=(2.0-c)*(x(n)-x(n1)) s(n)=s(n) - c*s(n1) c Back substitiute. s(n)= s(n)/a(n) do 4000 j=1,n1 i=n-j s(i) =(s(i) - (x(i+1)-x(i))*s(i+1))/a(i) 4000 continue return end c_____________________________________________________________________ function eval(y, nn, x, u, s) c$$$$ calls nothing c Performs cubic spline interpolation of a function sampled unequally c in x. The routine spline must be called to set up the array s. c See comments in spline. c y the coordinate at which function value is desired. c nn number of samples of original function. c x() nn-vector of x values in original data sample. c u() nn-vector of function values. c s() nn-vector of 2nd derivatives at sample points. Found by spline c which must be called once before beginning interpolation. c If y falls outside (x(1), x(nn)) extrapolation is based upon c the tangent straight line at the endpoint. common /startx/ istart dimension x(*),u(*),s(*) ce n=iabs(nn) if (istart.gt.n .or. istart.lt.1) $ istart=1+(n-1)*(y-x(1))/(x(n)-x(1)) if (y .le. x(1)) go to 3000 if (y .ge. x(n)) go to 3100 c Locate interval (x(l0),x(l1)) of y. if (y .lt. x(istart)) goto 1200 c Scan up the x array. do 1100 l1=istart, n if (x(l1) .gt. y) go to 1150 1100 continue 1150 l0=l1 - 1 go to 1500 c Scan downwards in x array. 1200 do 1300 l1=1, istart l0=istart - l1 if (x(l0) .le. y) go to 1350 1300 continue 1350 l1=l0 + 1 1500 istart=l0 c Evaluate interpolation. xi=(y - x(l0))/(x(l1) - x(l0)) h6=(x(l1) - x(l0))**2/6.0 eval=u(l0) + xi*(u(l1)-u(l0) - h6*(1.0-xi)*(s(l1)*(1.0+xi) $ + s(l0)*(2.0-xi))) return c Out of range. Substitute straight-line extrapolant. 3000 h=x(2) - x(1) eval=u(1) + (y-x(1))*((u(2)-u(1))/h - h*(s(2)-2.0*s(1))/6.0) return 3100 h=x(n) - x(n-1) eval=u(n)+(y-x(n))*((u(n)-u(n-1))/h+h*(2.0*s(n)-s(n-1))/6.0) return end c_____________________________________________________________________