clear all
close all
fclose all
hflnm='d:\1cva\bering\br';
[fnm,pthnm]=uigetfile('*.mat','Choose the IR .mat file');

a=['load ' pthnm fnm]
eval(a);

df=fr/ad;
[fnm,pnm]=uigetfile('*.grd','Choose the grid data file');
a=[pnm fnm];
[fgrd,msg]=fopen([pnm fnm],'rt');
if fgrd==-1,
   a,
   msg,
   return;
end
a=fgetl(fgrd)
nz=fscanf(fgrd,'%d',1)
ny=fscanf(fgrd,'%d',1)
%a=fgetl(fgrd)
zmin=fscanf(fgrd,'%g',1)
zmax=fscanf(fgrd,'%g',1)
dz=(zmax-zmin)/(nz-1);
ymin=fscanf(fgrd,'%g',1)
ymax=fscanf(fgrd,'%g',1)
dy=(ymax-ymin)/(ny-1);
cmin=fscanf(fgrd,'%g',1)
cmax=fscanf(fgrd,'%g',1)
a=fgetl(fgrd)
%[cwt,cnt]=fread(fgrd,'%g',[nz,ny]);
cwt=zeros(nz,ny);
for j=1:ny
   for i=1:nz
      cwt(i,j)=fscanf(fgrd,'%e ',1);
   end;
end;
fclose(fgrd);

[fnm,pnm]=uigetfile('*.dat','Choose the relief data file');
a=[pnm fnm];
[fgrd,msg]=fopen([pnm fnm],'rt');
if fgrd==-1,
   a,
   msg,
   return;
end
nxr=fscanf(fgrd,'%d',1)
a=fgetl(fgrd)
for i=1:nxr
   xrl(i)=fscanf(fgrd,'%e ',1);
   yrl(i)=fscanf(fgrd,'%e ',1);
end;
yrl(nxr+1)=ymax;xrl(nxr+1)=xrl(nxr);
yrl(nxr+2)=ymax;xrl(nxr+2)=zmax;
yrl(nxr+3)=ymin;xrl(nxr+3)=zmax;
fclose(fgrd);

fg5=figure;
set(gcf,'Units','centimeters');
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[1.0 1.0 9.6 11.0]);
subplot(411);
[c,h]=contour(ymin:dy:ymax,zmin:dz:zmax,cwt,1:max(max(cwt)),'k');
set(gca,'FontSize',10);
%title(pnm(size(pnm,2)-4:size(pnm,2)-1));
axis ij;
%xlabel('Range, km');
ylabel('Depth, m');
hold on;
line(yrl,xrl);
fill(yrl,xrl,'k');
clabel(c,h,'manual');
hold off;

subplot(412);
t=1/df;
dt=t/(2*(ad+1));
i=fix(1.0/dt);f=fix(2.0/dt);
shft=xr(1)/cgrmx
%contour(0:dz:dz*(nfn-1),i*dt:dt:f*dt,abs(ir(i+1:(f+1),:)))
contour((shft+i*dt):dt:(shft+f*dt),0:dz:dz*(nfn-1),abs(ir(i+1:(f+1),:))','k');
v=axis;
v=[20.6 21.8 v(3) v(4)];v=axis(v);
set(gca,'FontSize',10);
axis ij;
%xlabel('Time, sec');
ylabel('Depth, m');
%title('Complex envelope of the IR');
l=2*ad+1;
subplot(413);
plot((shft+i*dt):dt:(shft+f*dt),abs(ir(i+1:(f+1),22)),'k-');
v=axis;
%v=[0.4 1.0 0.0 30.0];v=[1.0 2.0 0.0 30.0];axis(v);
set(gca,'FontSize',10);
%xlabel('Time, sec');
ylabel('Rcv at 21 m');
%title('Complex envelope of the Impulse Response for 12 m depth receiver');
subplot(414);
plot((shft+i*dt):dt:(shft+f*dt),abs(ir(i+1:(f+1),41)),'k-');
v=axis;
%v=[1.0 2.0 0.0 30.0];axis(v);
set(gca,'FontSize',10);xlabel(['Time, sec                            ' pnm(size(pnm,2)-4:size(pnm,2)-1)]);
ylabel('Rcv at 40 m');
%title('Complex envelope of the Impulse Response for 23 m depth receiver');
zoom on;
a=['print -dtiff -loose ' [pnm 'fig5']]
eval(a);
%eval(a);
%disp('FFT done. Press any key...');pause;
if 0==1,
   fg1=figure;
      surf(0:dz:dz*(nfn-1),df*ifrlo:df:df*ifrhi,abs(pr))
   title(...
      'Amplitude of the transfer function versus frequency and depth');
   xlabel('Depth, km');ylabel('Frequency, Hz');
   whitebg;
   
   fg2=figure;
   subplot(211);
   plot(...
      df*ifrlo:df:df*ifrhi,abs(pr(:,23)),'b-');
   %df*ifrlo:df:df*ifrhi,abs(pr(:,23)),'gx',...
   xlabel('Frequency, Hz');ylabel('Amplitude');
   title(...
      'Transfer function for 22 m receiver depth');
   subplot(212);
   plot(...
      df*ifrlo:df:df*ifrhi,unwrap(angle(pr(:,23))),'b-');
   %df*ifrlo:df:df*ifrhi,unwrap(angle(pr(:,23))),'gx',...
   xlabel('Frequency, Hz');ylabel('Phase, grade');
   zoom on;
   %disp('Transfer function ready. Press any key...');pause;
   
   fg3=figure;
   t=1/df;
   dt=t/(2*(ad+1));
   i=0;f=2*ad+1;
   surf(0:dz:dz*(nfn-1),i*dt:dt:f*dt,abs(ir(i+1:(f+1),:)));
   xlabel('Depth, km');ylabel('Time, sec');
   title(...
      ['Amplitude of the complex envelope of the Impulse'...
         ' Response versus time and depth']);
   fg4=figure;
   t=1/df;
   dt=t/(2*(ad+1));
   i=0;f=1600;
   surf(0:dz:dz*(nfn-1),i*dt:dt:f*dt,abs(ir(i+1:(f+1),:)))
   xlabel('Depth, km');ylabel('Time, sec');
   title(...
      ['Amplitude of the complex envelope of the IR'...
         pthnm]);
end;

return

fg6=figure;
l=2*ad+1;
subplot(211);
plot(...
   0:dt:dt*l,abs(ir(:,13)),'b-')
%0:dt:dt*l,abs(ir(:,23)),'gx',...
xlabel('Time, sec');ylabel('Amplitude');
title(...
   'Complex envelope of the Impulse Response for 12 m depth receiver');
subplot(212);
%ph=unwrap(angle(ir(:,23).*exp(I*2*pi*400*dt*(0:l))'));
ph=unwrap(angle(ir(:,13)));
p=polyfit((0:dt:dt*l)',ph,1);
ph=ph-p(1)*(0:dt:dt*l)';
plot(...
   0:dt:dt*l,ph,'b-')
%0:dt:dt*l,ph,'gx',...
xlabel('Time, sec');ylabel('Phase, grade');
zoom on