% . . Test to see what angles should be plotted
clear all
close all
v=2000;
c=1120;
eps2=0.00001;

atodo=[0.1:0.1:179.9];
ttodo=[2:0.5:3.5]+0.04;

dz=10;
dx=10;
da=1.76471;
zaxis=[0:dz:5990];
xaxis=[-5000:dx:4990];
aaxis=[-45:da:43.2355];

ImageVolume=zeros(length(zaxis),length(xaxis),length(aaxis) );

for tt=1:length(ttodo)

  for ii=1:length(atodo)
    t = ttodo(tt);
    ang       = atodo(ii)/180*pi;
    ang1 (ii) = ang;
    x1        = 2/(1-4*c^2/(v^2*t^2));
    x2        = 1/(1-2*c*cos(ang)/(v*t));
    arg  (ii) = v*t/(2*c)*( 1 / (x1-x2) - 1);
    ang2 (ii) = acos(arg(ii));
    angle(ii) = (ang2(ii)-ang1(ii))/2;

    x(ii)=v*t/2*(cos(ang1(ii))+cos(ang2(ii)))/2;
    z(ii)=v*t/2*sqrt(1-c^2/v^2/t^2)*(sin(ang1(ii))+sin(ang2(ii)))/2;
    
    [ix,jx]=min(abs(x(ii)    -xaxis));
    [iz,jz]=min(abs(z(ii)    -zaxis));
    [ia,ja]=min(abs(-angle(ii)*180/pi-aaxis));
    ImageVolume(jz,jx,ja) = 1.;
  end
  
  figure(2)
  %ang2 = interp1(x,angle,[-2000:2000],'spline');
  hold on
  scatter(x,z,50,-angle*180/pi,'filled')
  axis ij  
  colorbar
  axis equal 
  axis tight
asdf
  
end

ta=xlabel('Distance (m)');
set(ta,'FontSize',12);
tb=ylabel('Depth (m)');
set(tb,'FontSize',12);

t2=text(5000,2800,'Reflection Angle (\circ)');
set(t2,'FontSize',12,'Rotation',90);

t1=text(-1120,-40,'\nabla');
set(t1,'FontSize',14);
t4=text(-1140,-200,'S');
set(t4,'FontSize',14);

t3=text( 1120,-40,'\nabla');
set(t3,'FontSize',14);
t5=text( 1100,-200,'R');
set(t5,'FontSize',14);

print -depsc2 Theory.ps

eval(['D1.history_file=''Theory.HH'' '])
D1.n=[length(zaxis) length(xaxis) length(aaxis)];
D1.o=[zaxis(1) xaxis(1) aaxis(1) ];
D1.d=[dz dx da];
D1.esize=4;
eval(['D1.in=''/net/pompei/scr1/jeff/Theory.HH@'' ']);
D1.data_format='xdr_float';
D1.data=ImageVolume;

eval(['sepwrite(D1,''Theory.HH'') '])
