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
gline=[30:0.25:150]/180*pi;
xtop=100*cos(gline)+9000; %  trayin(1:2:end);
ztop=100*sin(gline); %   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;

for ig=1:length(Rline)
  disp(['Working on Ray ',num2str(ig),' of ',num2str(length(Rline))])
  xx(1,ig)=xtop(ig);
  zz(1,ig)=ztop(ig);
  for ii=2:length(Rline)
    [nu,ix]=min(abs(xx(ii-1,ig) - xpts));
    [nu,iz]=min(abs(zz(ii-1,ig) - zpts));
    xx(ii,ig)=xx(ii-1,ig)+ds*FX(iz,ix)/K(iz,ix);
    zz(ii,ig)=zz(ii-1,ig)+ds*FZ(iz,ix)/K(iz,ix);
  end
  
end

etime1=TTimes1(iz,ix);

meshy=cat(3,xx,zz);
clear D1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% . . S velocity
% . . Create D1 Data file
   eval(['D1.history_file=''PS.rays.HH'' '])
   D1.n=[size(xx,1) size(xx,2) 2];   D1.o=[0.0 0.0 0];
   D1.d=[10 10 1];    D1.esize=4;
   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'') '])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
