%   |------------------------------------------------------------|
%   | $Id: delaytmR.m,v 1.2 2017/05/10 21:29:31 pm Exp $         |
%   | Basic Seismic Utilities (BSU)                              |
%   | LANGUAGE: octave                                           |
%   | AUTHOR: P. MICHAELS                                        |
%   | FUNCTION NAME: direct.m                                    |
%   | DIRECTORY: BSUROOT/octave                                  |
%   |  revised for OCTAVE                                        |
%   |------------------------------------------------------------|
%   | purpose: to compute delay time solution                    |
%   |   (normal refraction shooting, common shot gathers)        |
%   |Typical Flow: (pic first breaks, insert to headers)         |
%   | file.seg --> segpic.m -->pics.txt --| bpic                 |
%   |                                     |-->bpic-->bpicfile.seg|
%   |                          file.seg --|                      |
%   |                                                            |
%   | bpicfile.seg -->bref ---> G001 |                           |
%   |                      |--> D001 |-->direct.m --> overburden |
%   |                      |--> E001 |                velocity   |
%   |                                                & plot      |
%   | bpicfile.seg -->bref ---> G001 |                           |
%   |                      |--> D001 |-->delaytmR.m -->refractor |
%   |                      |--> E001 |                velocity   |
%   |                                                 structure  |
%   |------------------------------------------------------------|
%   | octave is an interpreter useful for scientifc calculations.|
%   | It is similar to MATLAB and may be obtained free from      |
%   | http://www.gnu.org/software/octave                         |
%   |------------------------------------------------------------|
%   |  NOTICE:                                                   |
%   |  Copyright (c) 2017 Paul Michaels                          |
%   |  <pm@cgiss.boisestate.edu>                                 |
%   |  This program is free software; you can                    |
%   |  redistribute it and/or modify it under the terms          |
%   |  of the GNU General Public License as published            |
%   |  by the Free Software Foundation; either version           |
%   |  3 of the License, or (at your option) any later           |
%   |  version.  This program is distributed in the              |
%   |  hope that it will be useful, but                          |
%   |  WITHOUT ANY WARRANTY; without even the implied            |
%   |  warranty of MERCHANTABILITY or FITNESS FOR A              |
%   |  PARTICULAR PURPOSE.  See the GNU General Public           |
%   |  License for more details.                                 |
%   |  You should have received a copy of the GNU                |
%   |  General Public License along with this program;           |
%   |  if not, write to the Free Software Foundation,            |
%   |  Inc., 675 Mass Ave, Cambridge, MA 02139, USA.             |
%   |------------------------------------------------------------|
clear
%fflush(stdout)
prompt = {'System Matrix File','Data Vector File','Elevations File',...
'Water Depth File'};
dlg_title = 'Enter bref files ';
num_lines = 1;
defaultans = {'G001','D001','E001','wds.dat'};
cn = inputdlg(prompt,dlg_title,num_lines,defaultans);
%frect = str2double(cn);
Gfile=cn{1,1}; Dfile=cn{2,1}; Efile=cn{3,1}; wdsfile=cn{4,1};

% read data
G=load('-ascii',Gfile);
D=load('-ascii',Dfile);
E=load('-ascii',Efile);
WDS=load('-ascii',wdsfile);

% assign variables
d=D(:,1);
flag=D(:,2);
elv=E(:,3);
chan=E(:,1);
recv=E(:,2);
wds=WDS(:,1);
iwds=WDS(:,2);

[N,M]=size(G);

prompt = {'Number of Receivers on Shore'};
dlg_title ='Enter Parameter';
num_lines =1;
defaultans ={num2str(1)};
cn=inputdlg(prompt,dlg_title,num_lines,defaultans);
frect=str2double(cn);
nshot=frect(1);  
nchanl=M-nshot-1;
pi2=pi/2;

prompt = {'Smoothness Weight 0<W<1'};
dlg_title = 'Enter Weight ';
num_lines = 1;
defaultans = {num2str(.1)};
cn = inputdlg(prompt,dlg_title,num_lines,defaultans);
frect = str2double(cn);
wsmooth=frect(1);
wsmtxt=sprintf('%.3f',wsmooth);
wfit=1-wsmooth;
for j=(nshot+1):(M-1)
for k=(nshot+1):M
W(j,k)=0;
end
W(j,j+1)=1;
W(j,j)=-1;
end
W(M,M)=1;
for j=1:(M-1)
W(M,j)=0;
end
%  Compute Smooth Inverse
H=inv(wsmooth*W'*W+wfit*G'*G)*wfit*G';

% Compute Delay Times
m=H*d;

for i=1:nshot
dltsht(i)=m(i);
end
k=1;
for i=(nshot+1):(M-1)
dltrcv(k)=m(i);
k=k+1;
end
v2=1/(m(M));
v2txt=sprintf('Refractor Velocity =%d m/s \n',ceil(v2));
sdtxt=sprintf('Receiver Delay Times (msec) \n');
sdtmtxt=sprintf('%.1f  ',dltsht*1000);
msgbox(strcat(v2txt,sdtxt,sdtmtxt));

% error analysis
cdd=eye(N,N);
sigmat=.001;  %ASSUME pics +/- uncertainty, 1 std dev (68% conf)
for j=1:N
cdd(j,j)=cdd(j,j)*sigmat*sigmat;
end
cm=H*cdd*H';
sigmaslow=sqrt(cm(M,M));
sigmav2=(1/(m(M)-sigmaslow)-1/(m(M)+sigmaslow) )/2;
k=0;
for j=(nshot+1):(M-1)
k=k+1;
errdt(k)=sqrt(cm(j,j));
end

%  plot delay times
figure(1)
clf
jrecv(1:length(dltrcv))=recv(1:length(dltrcv));
plot(jrecv,dltrcv*1000);
xmax=max(jrecv); xmin=min(jrecv);
ymax=max(dltrcv); ymin=min(dltrcv);

% label x-axis
xlabel('Shot Station Number');
% label y-axis
ylabel('Delay Time (msec)');
T=sprintf('Shot Delay Times (V2=%.0f m/s)',ceil(v2));
title(T)
grid on
disp('writing delay times to disk .  .  .');
fp2=fopen('dtimes.dat','w');
for j=1:length(dltrcv)
fprintf(fp2,'  %10.2f  %10.4f\n',jrecv(j),dltrcv(j));
end
fclose(fp2);
drawnow


prompt = {'Overburden Velocity (m/s)'};
dlg_title ='Enter Parameter';
num_lines =1;
defaultans ={num2str(500)};
cn=inputdlg(prompt,dlg_title,num_lines,defaultans);
frect=str2double(cn);
v1=frect(1);
  
% plot structural interpretation (overburden velocity constant)
theta=asin(v1/v2);
cs=cos(theta);
  if v2>1500 
  theta1=asin(1500/v2);  % p-refraction problem, ray optics
  else
  theta1=0.             % water over sv conversion problem, no ray optics
  end
  cst1=cos(theta1);     % cosine of angle in water layer
% h=(v1/cs)*dltrcv;
% calculate thickness of v1 layer, h
for jsp=1:length(dltrcv)
  if wds(jsp)>0. 
  h(jsp)=(v1/cs)*(dltrcv(jsp)-(wds(jsp)/cst1)/1500.);
  else
  h(jsp)=(v1/cs)*dltrcv(jsp);
  end
end % next jsp

NN=length(h);
for j=1:NN
e(j)=elv(j);
  if wds(j)>0. 
  ewds(j)=e(j)-wds(j);  % elevation of river bottom
  else
  ewds(j)=e(j);
  end
end
wx=ewds-h;

figure(2)
clf
%  plot structural solution
xmax=max(jrecv); xmin=min(jrecv);
ymax=max(e); ymin=min(wx);

plot(jrecv,wx,'-r');
hold on
plot(jrecv,e,'-b');
plot(iwds,ewds,'-b');   %use iwds to avoid order errors
grid on
axis([xmin,xmax,ymin,ymax]);
% label x-axis
xlabel('Shot Station Number');
% label y-axis
ylabel('Elevation (m)');
T=sprintf('Shot Structure (V1=%.0f   V2=%.0f m/s)',v1,v2);
title(T)

disp('writing structure to disk .  .  .');
fp2=fopen('struct.dat','w');
for j=1:length(e)
fprintf(fp2,'  %10.2f  %10.2f  %10.2f  %10.2f \n',...
jrecv(j),wx(j),e(j),ewds(j));
end
fclose(fp2)
drawnow

prompt = {'Alternative Fixed Depth Refractor (m)'};
dlg_title ='Enter Parameter';
num_lines =1;
defaultans ={num2str(10)};
cn=inputdlg(prompt,dlg_title,num_lines,defaultans);
frect=str2double(cn);
h1=frect(1);

figure(3)
clf
% plot velocity interpretation (refractor depth=constant)
for jr=1:NN
vovb(jr)=v2/sqrt(1+ (dltrcv(jr)*(v2/h1)).^2);
end
%  plot velocity solution
xmax=max(jrecv); xmin=min(jrecv);
ymax=max(vovb); ymin=min(vovb);

plot(jrecv,vovb,'-r');
h1txt=sprintf('V1 Velocity Solution (m/s) H1=%d m',ceil(h1));
title(h1txt);
% label x-axis
xlabel('Shot Station Number');
% label y-axis
ylabel('Overburden Velocity (m/s)');
axis([xmin,xmax,ymin,ymax]);
grid on

disp('writing velocity structure to disk   .  .');
fp2=fopen('velocity.dat','w');
for j=1:length(vovb)
fprintf(fp2,'  %10.2f  %10.2f \n',jrecv(j),vovb(j));
end
fclose(fp2);

% calculated times
dcal=G*m;

% least square error
e=d-dcal;
emin=sqrt(e'*e/length(e));
drawnow

% plot time fit
prompt = {'Number Constraint Equ. (if any)'};
dlg_title ='Enter Parameter';
num_lines =1;
defaultans ={num2str(0)};
cn=inputdlg(prompt,dlg_title,num_lines,defaultans);
frect=str2double(cn);
nconst=frect(1);
N2=(N-nconst);
% find break points between shots
jshot=1;
for j=1:N2
  if G(j,jshot)== 1 
  kend(jshot)=j;
  else
  jshot=jshot+1;
  end
end

figure(4)
clf
[nr,nc]=size(d);
nr2=nr-nconst; % index of last time, before constraints
xmin=min(flag);
xmax=max(flag);
ymin=min(min(d(1:nr2)),min(dcal(1:nr2)));
ymax=max(max(d(1:nr2)),max(dcal(1:nr2)));
kstart=1;
for k=1:nshot
if k==1 
  plot(flag(kstart:kend(k)),d(kstart:kend(k))*1000,'+');
hold on
else;
  plot(flag(kstart:kend(k)),d(kstart:kend(k))*1000,'+');
hold on
end;  % endif
  plot(flag(kstart:kend(k)),dcal(kstart:kend(k))*1000,'-r');
kstart=kend(k)+1;
end;  % next k
grid on
xlabel('STATION');
ylabel('Arrival Time (msec)');
axis([xmin,xmax,ymin*1000,ymax*1000]);
T=sprintf('+=Observed Times  --Solution')
title(T)

disp('writing data times to disk   .  .');
fp2=fopen('ttimes.dat','w');
fprintf(fp2,'  %5.5d  %5.5d \n',ceil(v1),ceil(v2));
fprintf(fp2,'  %5.5d \n',nshot);
for k=1:nshot
fprintf(fp2,' %5.5d \n',kend(k));
end
kendd=length(flag)-nconst;
for k=1:kendd
fprintf(fp2,'  %10.2f  %10.4f  %10.4f \n',flag(k),d(k),dcal(k));
end
fclose(fp2);

% write figures to postscript files  ONLY if in Octave, Not Matlab
%figure(1)
%print('delaytmplot.ps','-dpsc2');
%figure(2)
%print('structureplot.ps','-dpsc2');
%figure(3)
%print('velocityplot.ps','-dpsc2');
%figure(4)
%print('arrivaltimeplot.ps','-dpsc2');


