ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c taper.f is a toy Fortran program to play with tapers c c (a) gate function (no taper) c (b) triangular function (taper-1) c (b) Hanning function (taper-2) c c used in place of formal lecture notes p. 29-39 c started class with model prediction of c Arctic Ocean Surface Currents Jan.-1, 2016 c and more specifically near East Greenland c c Sep.-19, 2019 Andreas Muenchow, University of Delaware c compiles on OS X 10.14.6 with c GNU Fortran (GCC) 5.4.0 c c Oct.-06, 2021 Andreas Muenchow, University of Delaware c compiles on maxos 11.6 with c GNU Fortran (GCC) 11.2.0 c c runs in /Users/a2/data/copernicus c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc parameter (nn=155) real x(nn),y(nn),z(nn),z2(nn) real xr(nn),xi(nn),xr2,xi2 real arg,pi,freq(nn),Xspec real t,h double precision zvar,zvar2,Var integer n,k c pi = 4.*atan2(1.,1.) c c "ds" is 12.5 km grid spacing for "nn" pixel of the c European Commission's COPERNICUS c Marine Environmental Monitoring Service c Arctic Reanalysis Project at c http://marine.copernicus.eu c ds = 12.5 T0 = nn*ds c c Read Input Data and estimate Variance in Space Domain c open (unit=10,file='spd.tem') zvar = 0. do 1 i=1,nn read(10,*)x(i),y(i),z(i) c c Known Signal: 50 km oscillation c arg = 2.*pi/50.*float(i)*ds z(i) = 10.*cos(arg) c z(i) = 1. zvar = zvar+z(i)*z(i) 1 continue zvar = zvar/float(nn) write(6,*)'Variance of Input: ',zvar c c Apply Taper to minimize the sinc-function response of data Window c zvar2 = 0. do 20 i=1,nn t = -ds*float(nn/2)+float(i)*ds c c "z2" is the Do-Nothing-Taper (gate function) c h = 1. z2(i) = z(i)*h c c Triangular Taper: c c h = 1.-2.*sqrt(t*t)/T0 C z(i) = z(i)*h c z(i) = z(i)*h/sqrt(0.333442) c c Hanning Taper: c h = 0.5*(1.+cos(2.*pi*t/T0)) c z(i) = z(i)*h z(i) = z(i)*h/sqrt(0.375) zvar2 = zvar2+z(i)*z(i) c c write location "t", taper "h", data "z2" and tapered data "z: c write(12,*)t,h,z2(i),z(i) 20 continue zvar2 = zvar2/float(nn) write(6,*)'Variance of tapered input: ',zvar2 c cccccccccccccccccccccccccccccccccccccccccccccccccccc| c | c Fourier Transform Start | c V Var = 0. c c Loop over all frequencies "n" from 1 to nn/2 c (why nn/2 ? ---> Review Sampling Theorem) c do 2 n=1,nn/2 xr(n) = 0. xi(n) = 0. xr2 = 0. xi2 = 0. c c Loop over all data "k" from 1 to nn c do 3 k=1,nn c c "arg" is the argument of the exponential exp[-2*pi*k*n/nn] c "xr" and "xi" are real and imaginary parts of Fourier Transform c arg = -2.*pi*float(k*n)/float(nn) c tapered data xr(n) = xr(n) + z(k)*cos(arg) xi(n) = xi(n) + z(k)*sin(arg) c untapered data xr2 = xr2 + z2(k)*cos(arg) xi2 = xi2 + z2(k)*sin(arg) 3 continue c xr(n) = ds*xr(n) xi(n) = ds*xi(n) c xr2 = ds*xr2 xi2 = ds*xi2 c freq(n) = float(n)/T0 Xspec = 2.*(xr(n)*xr(n)+xi(n)*xi(n))/T0 Xspec2 = 2.*(xr2*xr2 +xi2*xi2) /T0 var = var+Xspec/T0 write(11,*)freq(n),Xspec,Xspec2 c 2 continue c c write Fourier Transform (complex) as c c negative frequency for imaginary part (sine) c write(13,101)(-freq(n),xi(n),n=nn/2,1,-1) c c positive frequency for real part (cosine) c write(13,101)(+freq(n),xr(n),n=1,nn/2,+1) c ^ c Fourier Transform End | c | cccccccccccccccccccccccccccccccccccccccccccccccccccc| Var = Var write(6,*)'Variance of Spectra: ',Var write(6,*)'Variance frequency/time: ',Var/zvar,Var/zvar2 101 format(2f20.6) stop end