
clear
close all

% This program (plus the subprogram "taper.m") reproduces some of the
% figures from  Igel, H., and O. Gudmundsson (1997), Frequency-dependent effects 
% on travel times and waveforms of long-period S waves: implications for the scale 
% of mantle heterogeneity, Phys. Earth Plan. Int., 104(1-3), 229-246.

clear

% dim
n=1024;
n2=1024;
ntap=40;
slope=input(' slope (e.g., -1, -2) ? ')
corner=input(' corner (e.g., 10, 20) ? ')


arff=(1:n)'*0.;
afff=(1:n/2)'*0.;

pause off


% taper in theta
h=1:(n2-2*ntap);
h=h*0+1;
h=taper(h',-50);
h=taper(h,50);
g=1:n2;
g=g*0;
g(ntap+1:n2-ntap)=h;


for i=1:20,

% initialize amplitude spectrum
q=(0:n/2-1);
q(1:corner)=(1:corner)*0+1;
for k=corner+1:n/2, q(k)=(k/(corner+.5))^(slope); end

% initialize phase
for k=1:n/2,
rF=q(k)*cos( 2 * pi * rand);                % Real part of spectrum
f(k)=rF+sqrt(-1)*sqrt(q(k)*q(k)-rF*rF);
end
for k=1:n/2-1
f(n+1-k)=conj(f(k+1));
end



ff=ifft(f);

rff=g.*real(ff);

%rff=taper(rff',-64);
%rff=taper(rff,64);


 subplot(311),loglog(q),title(' Idealized Spectrum '),pause

 subplot(312),plot(rff),title(' Perturbation '),pause

arff=arff+rff';

temp=abs(fft(rff));
temp=temp(1:n/2);


 subplot(313),loglog(temp),title(' Spectrum '),pause

afff=afff+temp';

end


afff=afff/max(afff);

m=100;
py=polyfit(log10(corner:m)',log10(afff(corner:m)),1);
nslope=py(1);
x=1:n/2;
xx=log10(corner:m);
yy=py(2)+py(1)*xx;

 subplot(311),loglog(q),title(sprintf(' Ideal spectrum: slope = %g',slope))
 subplot(312),plot(arff),title(' Summed Perturbation ')
 subplot(313),plot(log10(x),log10(afff),'-',xx,yy,'r-')
 title(sprintf(' Summed spectrum; nslope = %g', nslope))
