DRadiSurfTomo/utils/plot/plotgam.m
2025-01-14 12:02:34 +08:00

54 lines
1.2 KiB
Matlab

%mdlFile='DRadiSurfTomo.inMeasure.dat';
mdlFile='DRadiSurfTomo.inMeasure.dat.iter009';
plotDep=30; % plot horizantal slice at plotDep (30) km
dep=[0,5,10,15,20,25,30,35,40,45,50,55,60,70,80,90,100];
plotDepIndex=find(dep==plotDep);
ndep=length(dep);
minlat=28.0;
maxlat=32.50;
dlat=0.5;
lat=maxlat:-dlat:minlat;
nlat=length(lat);
minlon=90.5;
maxlon=98.5;
dlon=0.5;
lon=minlon:dlon:maxlon;
nlon=length(lon);
aa=load(mdlFile);
mdl=zeros(nlat,nlon,ndep);
% read the velocity model
i=1;
for idep=1:ndep
for ilon=1:nlon
for ilat=1:nlat
mdl(ilat,ilon,idep)=(aa(i,5)-1)/(aa(i,5)+1)*2;
i=i+1;
end
end
end
% smooth the velocity model
dlat=dlat/5;
dlon=dlat/5;
imagelat=maxlat:-dlat:minlat;
imagelon=minlon:dlon:maxlon;
[xin,yin]=meshgrid(lon,lat);
[xout,yout]=meshgrid(imagelon,imagelat);
image=griddata(xin,yin,squeeze(mdl(:,:,plotDepIndex)),xout,yout,'cubic');
% plot velocity
rd=[(0:31)/31,ones(1,32)];
gn=[(0:31)/31,(31:-1:0)/31];
bl=[ones(1,32),(31:-1:0)/31];
rwb=[rd',gn',bl'];
rwb=flipud(rwb);
imagesc(imagelon,imagelat,image); colormap(flipud(rwb)); colorbar('location','eastoutside');
%imagesc(imagelon,imagelat,image); colormap(flipud(jet)); colorbar('location','eastoutside');
hold on;
set(gca,'ydir','normal','Fontsize',14,'pos',[0.15 0.3 0.6 0.4]);