
pause on
clear
clf

ncol=64*3;
ncmp=3;cl=ncol/ncmp;
nx=201;
ny=201;
nz=201;
ns=4;
thres=.05;
thresbw=.3;
clip=.1;
iout=1;
trange=90;
prange=90;
t0=45;
p0=85;


dir

file=input(' Give name of run : ');



% Coordinates
disp(' Loading and initialising coordinates ')
eval(['load ' sprintf('%s_xco',file)]);theta=eval([sprintf('%s_xco',file)]);  
theta=theta(1:iout:nx);
eval(['load ' sprintf('%s_yco',file)]);phi=eval([sprintf('%s_yco',file)]);     
phi=phi(1:iout:ny);phi=phi-max(phi)/2;
eval(['load ' sprintf('%s_zco',file)]);r=eval([sprintf('%s_zco',file)]);  
r=r(1:iout:nz)/1000.;
nnx=length(theta);nny=length(phi);nnz=length(r);

% load PREM-Model
disp(' Loading and initialising model ')
load prem.dat; rp=6371-prem(:,1);vp=prem(:,3);
nvp=interp1(rp,vp,6371-r);
nvp(find(r>=6371))=min(vp);

nvp=nvp-min(nvp);nvp=nvp/max(nvp)*(cl-1)+1;
model=ones([1 nnz])'*nvp';


% load snaps
disp(' Which component ? ')
disp(' 1) P waves  ')
disp(' 2) shear (sides) and vertical (surface)')
disp(' 3) x component ')
disp(' 4) y component ')
disp(' 5) z component ')
ans=input(' Which ? ');


disp(' Loading data ... ')

if ans == 1,
eval(['load '  sprintf('%s_xvs',file) ]); snap1=eval([ sprintf('%s_xvs',file) ]); 
eval(['load '  sprintf('%s_yvs',file) ]); snap2=eval([ sprintf('%s_yvs',file) ]); 
eval(['load '  sprintf('%s_zvs',file) ]); snap3=eval([ sprintf('%s_zvs',file) ]); 
elseif ans == 2,
eval(['load '  sprintf('%s_xys',file) ]); snap1=eval([ sprintf('%s_xys',file) ]); 
eval(['load '  sprintf('%s_yxs',file) ]); snap2=eval([ sprintf('%s_yxs',file) ]); 
eval(['load '  sprintf('%s_zzs',file) ]); snap3=eval([ sprintf('%s_zzs',file) ]);    
elseif ans == 5,
eval(['load '  sprintf('%s_xzs',file) ]); snap1=eval([ sprintf('%s_xzs',file) ]); 
eval(['load '  sprintf('%s_yzs',file) ]); snap2=eval([ sprintf('%s_yzs',file) ]); 
eval(['load '  sprintf('%s_zzs',file) ]); snap3=eval([ sprintf('%s_zzs',file) ]); 
end
   
ns=length(snap1)/(nx*nz);

% amplification mode 
disp(' Which amplification ? ')
disp(' 0) None ')
disp(' 1) Sqrt ')
amp=input(' Which ? ');

% Threshold for snapshot display 
thres=input(' Threshold inside Earth ? ');

% Earth's surface

disp(' Initialising coordinates ')
r0=max(r);
x1=r0*sin(ones([1 nnx])'*theta').*cos(phi*ones([1 nny]));
y1=r0*sin(ones([1 nny])'*theta').*sin(phi*ones([1 nny]));
z1=r0*cos(ones([1 nnz])'*theta');

t=min(theta);
x2=ones([1 nnz])'*r'*sin(t).*cos(phi*ones([1 nny]));
y2=ones([1 nnz])'*r'*sin(t).*sin(phi*ones([1 nny]));
z2=ones([1 nnz])'*r'*cos(t);

p=min(phi);
x3=ones([1 nnz])'*r'.*sin(theta*ones([1 nny])).*cos(p);
y3=ones([1 nnz])'*r'.*sin(theta*ones([1 nny])).*sin(p);
z3=ones([1 nnz])'*r'.*cos(theta*ones([1 nny]));

t=pi/2;
x4=ones([1 nnz])'*r'*sin(t).*cos(phi*ones([1 nny]));
y4=ones([1 nnz])'*r'*sin(t).*sin(phi*ones([1 nny]));
z4=ones([1 nnz])'*r'*cos(t);

p=0;
x5=ones([1 nnz])'*r'.*sin(theta*ones([1 nny])).*cos(p);
y5=ones([1 nnz])'*r'.*sin(theta*ones([1 nny])).*sin(p);
z5=ones([1 nnz])'*r'.*cos(theta*ones([1 nny]));


% Earth's topography at the surface

disp(' Loading and initialising topography ')
load topo
[nt, np]=size(topo);
te=[0:nt-1]*180/nt;
pe=[0:np-1]*360/np;
nte=t0+[0:nnx-1]*trange/nnx;
npe=p0+[0:nny-1]*prange/nny;

topo=topo-min(min(topo));
topo=topo/max(max(topo))*cl+.55;
ntopo=interp2(pe,te',topo,npe,nte');


% initialise colormaps

disp(' Initialising colormaps ')
colormap(topomap1);
map1=colormap;

colormap('default');
colormap('gray');
map2=colormap;

colormap('default');
colormap(ry(.4));
map3=colormap;

map=[map1;map2;map3];


is0=input(' First snap ? ');
isinc=input(' Snap increment ? ');





for is=is0:isinc:ns;
   
   disp(sprintf(' Preparing snap %g ',is)); 
   
   % surface

f1=ntopo;
s=snap3( (is-1)*nz*ny+1:is*nz*ny );
s=reshape(s,nx,ny);
s=s(1:iout:ny,1:iout:nz)';
s=s/max(max(abs(s)));
s=enhance(s,thresbw,clip);
f1(find(s>thresbw)) = s(find(s>thresbw))*cl/2+1.5*cl;
f1(find(s<-thresbw))= s(find(s<-thresbw))*cl/2+1.5*cl;


% theta const

f2=model+cl;

% phi = cons

f3=model+cl;


% the two inner faces

% theta = const
f4=model+cl;

% the snap
s=snap1( (is-1)*nz*ny+1:is*nz*ny );
s=reshape(s,nz,ny);
s=s(1:iout:ny,1:iout:nz)';
s=s/max(max(abs(s)));

if amp == 1, s = (abs(s)).^(0.5).*sign(s); end

s=enhance(s,thres,clip);
f4(find(s>thres)) = s(find(s>thres))*cl/2+2.5*cl;
f4(find(s<-thres))= s(find(s<-thres))*cl/2+2.5*cl;

% phi = cons
f5=model+cl;
% the snap
s=snap2( (is-1)*nz*ny+1:is*nz*ny );
s=reshape(s,nz,ny);
s=s(1:iout:ny,1:iout:nz)';
s=s/max(max(abs(s)));

if amp == 1, s = (abs(s)).^(0.5).*sign(s); end

s=enhance(s,thres,clip);
f5(find(s>thres)) = s(find(s>thres))*cl/2+2.5*cl;
f5(find(s<-thres))= s(find(s<-thres))*cl/2+2.5*cl;

nnx2=round(nnx/2);
hold on
surf(-z1(nnx2:nnx,:),y1(nnx2:nnx,:),x1(nnx2:nnx,:),f1(nnx2:nnx,:))
surf(-z1(:,nnx2:nnx),y1(:,nnx2:nnx),x1(:,nnx2:nnx),f1(:,nnx2:nnx))
surf(-z2(nnx2:nnx,:),y2(nnx2:nnx,:),x2(nnx2:nnx,:),f2(nnx2:nnx,:))
surf(-z3(nnx2:nnx,:),y3(nnx2:nnx,:),x3(nnx2:nnx,:),f3(nnx2:nnx,:))
surf(-z4(1:nnx2+1,:),y4(1:nnx2+1,:),x4(1:nnx2+1,:),f4(1:nnx2+1,:))
surf(-z5(1:nnx2+1,:),y5(1:nnx2+1,:),x5(1:nnx2+1,:),f5(1:nnx2+1,:))
axis('equal'),
axis('off')

%view([-45 40])
view(3)
colormap(map);caxis([1 ncol]);
shading interp
if is==1,
light('Position',[-3 -1 3])
light('Position',[1 3 2])
light('Position',[-1 -1 -1])
end
pause
hold off

disp(sprintf(' This is snap number : %g',is))
ans=input(' Plot this snap ? ')
if ans == 1,
   fileout=input(' Give filename *.jpg ');
   eval(['print -djpeg100 ' fileout]);
end

end

