clear all ; close all ; clc

% . . Elliptic test plotting...
A = load_history('SHOTtest4.H');
B = load_data(A);
C = B.data;
C = C/abs(max(C(:)));

zaxis = ([0:B.n(1)]*B.d(1)+B.o(1))/1000.;
xaxis = ([0:B.n(2)]*B.d(2)+B.o(2))/1000.;
aaxis = [0:B.n(3)]*B.d(3)+B.o(3);

acut=20;

runIV=1;
if (runIV ==1 )
	v=2;
	c=1.120;
	eps2=0.00001;

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

	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
	end

	save ImageVolume.mat ImageVolume

else
    load ImageVolume.mat
end

colormap('gray')

clear A B C;

% . . For elliptic
A = load_history('SHOTtest3.H');
B = load_data(A);
C = B.data;

C = C/abs(max(C(:)));

axes('position',[0.075 0.56 0.5 0.4])
imagesc(xaxis,zaxis,C(:,:,acut))
set(gca,'YTickLabel',{'0.0'; '0.5'; '1.0'; '1.5'; '2.0'; '2.5'; '3.0'; '3.5'; '4.0'})
%l1=text(-3.5,-0.35,'a)');
%set(l1,'FontSize',12)

[ii,jj]=find(ImageVolume(:,:,acut));
for iii=1:length(ii)
	t1=text(xaxis(jj(iii)),zaxis(ii(iii)),'\bullet');
	set(t1,'FontSize',6);
end
clear ii jj;

t1=ylabel('Depth (km)');
t2=xlabel('Distance (km)');
set(t1,'FontSize',12)
set(t2,'FontSize',12)

c1=2/8;
c2=3/8;
c3=3/8;

t1=line([xaxis( floor(c1*B.n(2)))  xaxis( floor(c1*B.n(2)))],[0 4]);
set(t1,'LineWidth',2);
t1=line([ xaxis( floor(c2*B.n(2)))  xaxis( floor(c2*B.n(2)))],[0 4]);
set(t1,'LineWidth',2);
t1=line([xaxis( floor(c3*B.n(2)))  xaxis( floor(c3*B.n(2)))],[0 4]);
set(t1,'LineWidth',2);


axes('position',[0.635 0.56 0.08 0.4])
imagesc(aaxis,zaxis,squeeze(C(:,c1*B.n(2),:)) )
t1=ylabel('Depth (km)');
t2=xlabel('Angle (\circ)');
set(t1,'FontSize',10);
set(t2,'FontSize',10);
set(gca,'YTickLabel',{'0.0'; '0.5'; '1.0'; '1.5'; '2.0'; '2.5'; '3.0'; '3.5'; '4.0'})
%l1=text(-75,-0.35,'b)');
%set(l1,'FontSize',12)
[ii,jj]=find(ImageVolume(:,c1*B.n(2),:));
for iii=1:length(ii)
	t1=text(aaxis(jj(iii)),zaxis(ii(iii)),'\bullet');
	set(t1,'FontSize',10);
end
clear ii jj;


axes('position',[0.7775 0.56 0.08 0.4])
imagesc(aaxis,zaxis,squeeze(C(:,c2*B.n(2),:)) )
t1=ylabel('Depth (km)');
t2=xlabel('Angle (\circ)');
set(t1,'FontSize',10);
set(t2,'FontSize',10);
set(gca,'YTickLabel',{'0.0'; '0.5'; '1.0'; '1.5'; '2.0'; '2.5'; '3.0'; '3.5'; '4.0'})
%l1=text(-75,-0.35,'c)');
%set(l1,'FontSize',12)
[ii,jj]=find(ImageVolume(:,c2*B.n(2),:));
for iii=1:length(ii)
	t1=text(aaxis(jj(iii)),zaxis(ii(iii)),'\bullet');
	set(t1,'FontSize',10);
end
%clear ii jj;

%axes('position',[0.9175 0.56 0.08 0.4])
%imagesc(aaxis,zaxis,squeeze(C(:,c3*B.n(2),:)) )
%t1=ylabel('Depth (km)');
%t2=xlabel('Angle (\circ)');
%set(t1,'FontSize',10);
%set(t2,'FontSize',10);
%set(gca,'YTickLabel',{'0.0'; '0.5'; '1.0'; '1.5'; '2.0'; '2.5'; '3.0'; '3.5'; '4.0'})
%%l1=text(-75,-0.35,'d)');
%set(l1,'FontSize',12)
%[ii,jj]=find(ImageVolume(:,c3*B.n(2),:));
%for iii=1:length(ii)
%	t1=text(aaxis(jj(iii)),zaxis(ii(iii)),'\bullet');
%	set(t1,'FontSize',10);
%end

%brighten(0.5)
print -depsc2 ETtop.eps
disp(['Finished writing image file'])

