%
% get amplification of GL-1M/GL-100M and GL-1M/GL-40M at ETTYUUJIMA
% data is before 1985 to except the effects of SINDOUTOU
% 
% case 2: ( average of exception of amp )
%
% setup graph

fig3=figure('Units','centimeters','Position',[5.0 1.0 24.0 15.0]...
,'PaperUnits','centimeters','PaperPosition',[1.0 1.0 24.0 15.0]);
grh3=mysubplot(2,3,0,10,'Hz','cm/s/s*sec');
fig4=figure('Units','centimeters','Position',[5.0 1.0 24.0 15.0]...
,'PaperUnits','centimeters','PaperPosition',[1.0 1.0 24.0 15.0]);
grh4=mysubplot(2,3,0,10,'Hz','cm/s/s*sec');


% read data

fn='Macintosh_HD:Applications:Math:MATLAB 4.2c:toolbox:Mytool:ettyujima:jishin';
[obs,time,namm,dt,n,code]=wvrdf2(fn);

% calculate amplification

t=0:dt:dt*(n-1);
nyq=n/2+1;
df=1/(n*dt);
f=0:df:df*(nyq-1);
nwv=size(obs,2);
ndr=3;
ndp=3;
nev=round(nwv/ndr/ndp);

mgn=[ 1.0 1.0 1.0 1.0 ];

	spc=zeros(n,ndr*ndp);
	a100=zeros(nyq,ndr);
	a040=zeros(nyq,ndr);
	b100=zeros(nyq,ndr);
	b040=zeros(nyq,ndr);

for iev=1:nev

	fig1=figure('Units','centimeters','Position',[1.0 1.0 24.0 15.0]...
	,'PaperUnits','centimeters','PaperPosition',[1.0 1.0 24.0 15.0]);
	grh1=mysubplot(2,3,0,20,'sec','cm/s/s');
	is=ndr*ndp*(iev-1)+1;

% foruier transform

	for j=1:ndr
		for i=1:ndp

			ij=ndr*(i-1)+j;
			k=is+ij-1;
			fprintf('%20s \n',namm(k,1:20));
			spc(:,ij)=abs(fft(obs(:,k)));
			
		end;

			k100=is+j  -1;
			k040=is+3+j-1;
			k000=is+6+j-1;
			fprintf('%20s / %20s\n',namm(k000,1:20),namm(k040,1:20));
			a040(:,j)=a040(:,j)+spc(1:nyq,6+j)./spc(1:nyq,3+j);
			fprintf('%20s / %20s\n',namm(k000,1:20),namm(k100,1:20))
			a100(:,j)=a100(:,j)+spc(1:nyq,6+j)./spc(1:nyq,  j);

			figure(fig1);subplot(grh1(1,j));hold on;
			axis([0.0 10.0 0.0 5.0]);hold on;	
			set(gca,'title',text('String',namm(k040,1:20),'FontSize',8));
			plot(f,spc(1:nyq,6+j)./spc(1:nyq,3+j));
			
			figure(fig1);subplot(grh1(2,j));hold on;
			axis([0.0 10.0 0.0 5.0]);hold on;	
			set(gca,'title',text('String',namm(k100,1:20),'FontSize',8));
			plot(f,spc(1:nyq,6+j)./spc(1:nyq,  j));
				
	end;

%	pause;
	delete(fig1);
end;

for idr=1:ndr

	if idr == 1
		dir='NS';
	elseif idr == 2
		dir='EW';
	else
		dir='UD';
	end

	a040(:,idr)=a040(:,idr)/nev;
	a100(:,idr)=a100(:,idr)/nev;
	b040=parzn_etc(a040);
	b100=parzn_etc(a100);

		figure(fig3);
		subplot(grh3(1,idr));hold on;
		titl= sprintf('GL-1M/GL-40M %2s',dir);
		set(gca,'title',text('String',titl,'FontSize',8));
		plot(f,a040(:,idr),'r');
		
		subplot(grh3(2,idr));hold on;
		titl= sprintf('GL-1M/GL-100M %2s',dir);
		set(gca,'title',text('String',titl,'FontSize',8));
		plot(f,a100(:,idr),'r');
	
		figure(fig4);
		subplot(grh4(1,idr));hold on;
		titl= sprintf('SMOOTHED GL-1M/GL-40M %2s',dir);
		set(gca,'title',text('String',titl,'FontSize',8));
		plot(f,b040(:,idr),'w');
		
		subplot(grh4(2,idr));hold on;
		titl= sprintf('SMOOTHED GL-1M/GL-100M %2s',dir);
		set(gca,'title',text('String',titl,'FontSize',8));
		plot(f,b100(:,idr),'w');

		tt=sprintf('%2s',dir);
		titlk(1:2,idr)=tt';
end;

%output kaleida-graph

fn='Macintosh_HD:etc040.2.amp';
y = out_etc(fn,f,b040,titlk);

fn='Macintosh_HD:etc100.2.amp';
y = out_etc(fn,f,b100,titlk);

fprintf('WAVES : %5d \n',nev);

pause;

%delete(fig0);
%delete(fig2);
delete(fig3);
delete(fig4);
