clear all;
cd /net/koko/homes/sep/jeff/
startup;
cd /net/koko/homes/sep/jeff/Research/GEO_NO/
tmax=1.5;

file_name='PSeikonal.H';
D1=load_history(file_name);
fidd1=fopen(D1.in,'r','ieee-be');
arr_size=[D1.n(1),D1.n(2)];
nelem=arr_size(1)*arr_size(2);
f(1:nelem) = fread(fidd1,nelem,'float32');
TTimes=reshape(f,[D1.n(1),D1.n(2)]);

xpts=[0:D1.n(2)-1]*D1.d(2)+D1.o(2);
zpts=[0:D1.n(1)-1]*D1.d(1)+D1.o(1);

% . . Travel Time Tables
TTimes1=TTimes-min(min(TTimes));

D2=D1;
clear D1 f arr_size nelem 
file_name='rayin.H';
D1=load_history(file_name);
fidd1=fopen(D1.in,'r','ieee-be');
arr_size=2*D1.n;
nelem=arr_size(1);
trayin = fread(fidd1,nelem,'float32');

% . . Top Points
xtop=trayin(1:2:end);
ztop=trayin(2:2:end);

% . . For Side calculations
[FX,FZ]=gradient(TTimes1);
K=sqrt(FX.*FX+FZ.*FZ);
Rline=linspace(0.01,tmax,length(xtop));
ds=5;

% . . Left Points
xlft(1)=xtop(1);
zlft(1)=ztop(1);
for ii=2:length(Rline)
  [nu,ix]=min(abs(xlft(ii-1) - xpts));
  [nu,iz]=min(abs(zlft(ii-1) - zpts));
  xlft(ii)=xlft(ii-1)+ds*FX(iz,ix)/K(iz,ix);
  zlft(ii)=zlft(ii-1)+ds*FZ(iz,ix)/K(iz,ix);
end
etime1=TTimes1(iz,ix);

% . . Right Points
xrgt(1)=xtop(end);
zrgt(1)=ztop(end);
for ii=2:length(Rline)
  [nu,ix]=min(abs(xrgt(ii-1) - xpts));
  [nu,iz]=min(abs(zrgt(ii-1) - zpts));
  xrgt(ii)=xrgt(ii-1)+ds*FX(iz,ix)/K(iz,ix);
  zrgt(ii)=zrgt(ii-1)+ds*FZ(iz,ix)/K(iz,ix);
end
etime2=TTimes1(iz,ix);

ntmax=(etime1+etime2)/2;

% . . Bottom Points
[c,h]=contour(TTimes1,[ntmax ntmax]);
close all;
xr(1:length(c)-1) =(c(1,2:end)-1)*D2.d(2)+D2.o(2);
zr(1:length(c)-1) =(c(2,2:end)-1)*D2.d(1)+D2.o(1);

[nu,iil]=min( (xlft(end) - xr).^2+(zlft(end)-zr).^2) ;
[nu,iir]=min( (xrgt(end) - xr).^2+(zrgt(end)-zr).^2) ;
lall=length(iil:iir);

ccc(2:lall-1)=sqrt( (xr(iil+2:iir)-xr(iil:iir-2)).^2+(zr(iil+2:iir)- ...
                                                  zr(iil:iir-2)).^2)/2;
ccc(1)=ccc(2);
ccc(lall)=ccc(lall-1);

cccnorm=cumsum(ccc)/sum(ccc);

nxline=linspace(0,1,length(xtop));
Nxbot=interp1(cccnorm,xr(iil:iir),nxline,'spline');
Nzbot=interp1(cccnorm,zr(iil:iir),nxline,'spline');

figure(1)
plot(xtop,ztop,'r.',xrgt,zrgt,'b.',Nxbot,Nzbot,'k.',xlft,zlft,'g.')
axis ij

len=length(ztop);

x(1:len,1)=xlft;
z(1:len,1)=zlft;

x(1:len,len)=xrgt;
z(1:len,len)=zrgt;

x(1,1:len)=xtop;
z(1,1:len)=ztop;

x(len,1:len)=Nxbot;
z(len,1:len)=Nzbot;

meshy=complex(x,z);

dt=abs(mean(diff(xlft)));

clear D1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% . . S velocity
% . . Create D1 Data file
   eval(['D1.history_file=''PS.rays.HH'' '])
   D1.n=[len len 1];   D1.o=[0.0 0.0 0];
   D1.d=[dt dt 1];    D1.esize=8;
   eval(['D1.in=''/net/koko/jeff/PS.rays@'' ']);
   D1.data_format='xdr_float';
   D1.data=meshy;

% . . Call sepwrite program ;
   eval(['sepwrite(D1,''PS.rays.HH'') '])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Rangle=atan2(zbot,xbot-9000)/pi;
%nxline=linspace(min(Rangle),max(Rangle),length(Rangle));
%Nxbot=interp1(Rangle,xbot,fliplr(nxline));
%Nzbot=interp1(Rangle,zbot,fliplr(nxline));
%x2(1:len,1)=xlft;
%z2(1:len,1)=zlft;
%x2(1:len,len)=xrgt;
%z2(1:len,len)=zrgt;
%x2(1,1:len)=xtop;
%z2(1,1:len)=ztop;
%x2(len,2:len-1)=Nxbot(2:len-1);
%z2(len,2:len-1)=Nzbot(2:len-1);
%meshy2=complex(x2,z2);
%dt=abs(mean(diff(xlft)));
%clear D1
%eval(['D1.history_file=''Srays2.HH'' '])
%D1.n=[len len 1];   D1.o=[0.0 0.0 0];
%D1.d=[dt dt 1];    D1.esize=8;
%eval(['D1.in=''/net/koko/jeff/Srays2@'' ']);
%D1.data_format='xdr_float';
%D1.data=meshy2;
%eval(['sepwrite(D1,''Srays2.HH'') '])








