clear all; clc;
% Cartesian Test
nx=512;
dx=7.5;
ox=4000.;
x0=28000;

ny=1;
dy=20;
oy=0.;

nz=401;
dz=2.4*pi/(3*(nz-1));
oz=0.;

xi1= -fliplr([0:nx-1]*dx+ox);
xi2= [0:ny-1]*dy+oy;
xi3= ([0:nz-1]*dz+oz);

A=1.*ones(size(xi3)) +.3*xi3 -.05*xi3.^2;
%A=1.*ones(size(xi3))+( .4*sin(xi3)  )+ .2*xi3;

% Say v=v(z)  = 1500 + z
for ix=1:nx
  for iy=1:ny
    for iz=1:nz
      rays(iz,ix,iy,1)= A(iz)*xi1(ix)*cos(xi3(iz)) + x0; % X Coordinate
      rays(iz,ix,iy,2)= xi2(iy); % Y Coordinate
      rays(iz,ix,iy,3)= A(iz)*xi1(ix)*sin(xi3(iz)); % Z Coordinate
      Vel (iz,ix,iy)  = 1500 + 0.2*rays(iz,ix,iy,3);
    end
  end
end


outcoord=1;

if outcoord==0;

  pcut=8;
  figure(2); clf
  for iz=1:10:nz
    plot( squeeze(rays(iz,1:pcut:end,1,1)),...
         -squeeze(rays(iz,1:pcut:end,1,3)),'r');
    hold on
  end
  
  for ix=1:10:nx
    plot( squeeze(rays(1:pcut:nz,ix,1,1)),...
         -squeeze(rays(1:pcut:nz,ix,1,3)),'b');
    hold on
  end
  
  plot(squeeze(rays(:, 1, 1,1))',...
      -squeeze(rays(:, 1, 1,3))','r');
  plot(squeeze(rays(:,nx, 1,1))',...
      -squeeze(rays(:,nx, 1,3))','r');
  axis ij
  axis([15000 42000 0 12000])
end


  dz=mean(mean(mean(sqrt( ...
         ( rays(2:nz,:,:,1)-rays(1:nz-1,:,:,1) ).^2+...
         ( rays(2:nz,:,:,2)-rays(1:nz-1,:,:,2) ).^2+...
         ( rays(2:nz,:,:,3)-rays(1:nz-1,:,:,3) ).^2))));

if outcoord==1;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % . . Coordinate system
  eval(['D1.history_file=''Test1.rays.HH'' '])
  D1.n=[nz nx ny 3];   D1.o=[oz ox oy 0];
  D1.d=[dz dx dy 1];   D1.esize=4;
  eval(['D1.in=''/net/koko/jeff/Test1.rays@'' ']);
  D1.data_format='xdr_float';
  D1.data=rays;
  eval(['sepwrite(D1,''Test1.rays.HH'') '])
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 clear D1
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % . . Velocity
  eval(['D1.history_file=''Test1.Vel.HH'' '])
  D1.n=[nz nx ny];   D1.o=[oz ox oy];
  D1.d=[dz dx dy];   D1.esize=4;
  eval(['D1.in=''/net/koko/jeff/Test1.Vel@'' ']);
  D1.data_format='xdr_float';
  D1.data=Vel;
  eval(['sepwrite(D1,''Test1.Vel.HH'') '])
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
