clear all; close all; clc;
% Cartesian Test
nx=256;
dx=10;
ox=0.;

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

nz=150;
dz=10;
oz=0.;

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

theta=65/180*pi;

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


outcoord=1;

if outcoord==0;
  pcut=8;
  figure(1); clf
  for iz=1:10:nz
    plot3( squeeze(rays(iz,1:pcut:end,1:pcut:end,1)),...
           squeeze(rays(iz,1:pcut:end,1:pcut:end,2)),...
           -squeeze(rays(iz,1:pcut:end,1:pcut:end,3)),'r');
    hold on
    plot3(squeeze(rays(iz,1:pcut:end,1:pcut:end,1))',...
          squeeze(rays(iz,1:pcut:end,1:pcut:end,2))',...
         -squeeze(rays(iz,1:pcut:end,1:pcut:end,3))','b');
  end
  plot3(squeeze(rays(:,1,1,1))',...
        squeeze(rays(:,1,1,2))',...
        -squeeze(rays(:,1,1,3))','b');
  plot3(squeeze(rays(:,nx,1,1))',...
        squeeze(rays(:,nx,1,2))',...
        -squeeze(rays(:,nx,1,3))','b');
  plot3(squeeze(rays(:,1,ny,1))',...
        squeeze(rays(:,1,ny,2))',...
        -squeeze(rays(:,1,ny,3))','b');
  plot3(squeeze(rays(:,nx,ny,1))',...
        squeeze(rays(:,nx,ny,2))',...
        -squeeze(rays(:,nx,ny,3))','b');
  view([-4,8])
end
axis vis3d

  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=''Test0.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/Test0.rays@'' ']);
  D1.data_format='xdr_float';
  D1.data=rays;
  eval(['sepwrite(D1,''Test0.rays.HH'') '])
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 clear D1
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % . . Velocity
  eval(['D1.history_file=''Vel0.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/Vel0@'' ']);
  D1.data_format='xdr_float';
  D1.data=Vel;
  eval(['sepwrite(D1,''Vel0.HH'') '])
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
