clear all;
tmax=2.25;
R=10;
cd ~jeff/;startup;cd ~jeff/Research/GEO_NO;

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)]);

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

D2=D1;
%clear D1 f arr_size nelem 
%file_name='Geikon.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');

x=R*sin( ([-80:0.5:80])/180*pi)+9000;
z=R*cos( ([-80:0.5:80])/180*pi)+D2.o(1);

% . . Top Points
xtop=x(1:end);
ztop=z(1:end);

% . . Bottom Points
[c,h]=contour(TTimes1,[tmax tmax]);
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);

xline=[0:length(xtop)-1]/(length(xtop)-1);
xin  =[0:length(xr  )-1]/(length(xr  )-1);
xbot=interp1(xin,xr,xline);

zline= [0:length(ztop)-1]/(length(ztop)-1);
zin  = [0:length(zr  )-1]/(length(zr  )-1);
zbot = interp1(zin,zr,zline);


% . . Right Side Points
zrgt=zeros(1,length(zbot));
orx=xtop(end);
drx=(xbot(end)-xtop(end))/(length(xbot)-1);
xrgt=[0:length(xbot)-1]*drx+orx;

% . . Left Side Points
zlft=zeros(1,length(zbot));
olx=xtop(1);
dlx=(xbot(1)-xtop(1))/(length(xbot)-1);
xlft=[0:length(xbot)-1]*dlx+olx;

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

len=length(zbot);

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)=xbot;
z(len,1:len)=zbot;

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)/pi;

Nxbot=interp1(Rangle,xbot,fliplr(xline));
Nzbot=interp1(Rangle,zbot,fliplr(xline));

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=''PS.rays2.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.rays2@'' ']);
D1.data_format='xdr_float';
D1.data=meshy2;

eval(['sepwrite(D1,''PS.rays2.HH'') '])








