Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
143 lines (123 sloc) 4.17 KB
function varargout=readJRmodel(in,rmav,degres,toplot)
% [cosi,dels,dems]=readJRmodel
% lmcosi=readJRmodel(in,rmav)
% [r,lon,lat]=readJRmodel(in,rmav,degres)
% [r,lon,lat,ch,ph,cb,xl]=readJRmodel(in,rmav,degres,1)
% ...=readJRmodel(mname)
%
% Reads/converts/plots the spherical harmonic coefficients of the seismic
% wavespeed model S40RTS in the format by Jeroen Ritsema.
%
% INPUT: (if no input, returns the entire model over all splines)
%
% in The number of the spline "layer" function [default: 3]
% rmav Remove the degree-zero contribution [default: 0]
% degres Degree resolution of the expansion [default: 1]
% toplot Make a Mollweide map [default: 0]
% mname The model name {'S40RTS.sph' [default] | 's40_prelim.sph'}
%
% OUTPUT:
%
% cosi,dels, dems The coefficients for all of the splines/layers, OR
% lmcosi [degr order cos sin] coefficients for a single layer, OR
% r, lon,lat The map values and the longitudes and latitudes AND
% ch,ph,cb,xl Coninent, plate, colorbar, and xlabel handles of the figure
%
% SEE ALSO: PLM2XYZ, PLOTPLM, INTERPJRMODEL, RTSSPLINES, CALCJRMODEL
%
% EXAMPLE:
%
% readJRmodel(3,[],1,1) % should be the same as
% [a,dels,dems]=readJRmodel; plotplm([dels dems a(:,:,3)],[],[],[],1);
% kelicol % Note that this plots a single spline, which isn't yet a "depth"
%
% Last modified by fjsimons-at-alum.mit.edu, 11/10/2016
% Define default file names and directory locations
defval('diro',fullfile(getenv('IFILES'),'EARTHMODELS','RITSEMA'))
if exist('in','var') && ~isempty(in) && isstr(in)
% First input was a model name; this is a hack for calcJRmodel
mname=in;
% Sad hack for calcJRmodel which was replaced by interpJRmodel
flag=1;
end
defval('mname','S40RTS.sph')
% Open and load the file and read it in
fid=fopen(fullfile(diro,mname));
% Read the first line
fline=fgetl(fid);
% Will parse this line later but now I just "know" it
L=40; n=21;
% Read all data in one go - not sure why it doesn't like %11e4
cofs=reshape(fscanf(fid,'%f'),(L+1)^2,n);
fclose(fid);
% Check that the size is what we expect
difer(prod(size(cofs))-(L+1)^2*n,[],[],NaN)
% Create the vectors with harmonic degrees and orders, etc
[dems,dels,mz,lmcosi,mzi,mzo]=addmon(L);
% Stick the coefficients in at the right place
cosi=repmat(lmcosi(:,3:4),[1 1 n]);
% We're doing this in a roundabout way but the result counts
coso=cosi(:,:,1);
for ind=1:n
coso(mzo)=cofs(:,ind);
cosi(:,:,ind)=coso;
end
% Remove the spherical average of the RELATIVE PERTURBATION also
sphav=squeeze(cosi(1,1,:));
defval('rmav',0)
if rmav==1
cosi(1,1,:)=0;
end
% Need to include the (-1)^m factor also
CSC=(-1).^dems;
% A further adjustment that is needed and was thoroughly verified
dom=1./sqrt(2-(dems==0))/sqrt(4*pi);
if nargin>0 && flag~=1
% Identify the spline/layer to plot
defval('in',4);
% Collect the sine and cosine coefficients containing the model
plotcosi=cosi(:,:,in).*repmat(CSC.*dom,1,2);
% Change the relative to the PERCENTAGE perturbation
plotcosi=plotcosi*100;
% Prepare for output
varns={[dels dems plotcosi]};
elseif nargin==0 || flag==1
% if nothing else requested return all coefficients
for in=1:size(cosi,3)
cosi(:,:,in)=cosi(:,:,in).*repmat(CSC.*dom,1,2)*100;
end
% Prepare for output
varns={cosi,dels,dems};
end
% Output this if there are no further requests
if nargin>2 && flag~=1
% Perform the expansion
defval('degres',1)
[modl,lon,lat]=plm2xyz([dels dems plotcosi],degres);
% Prepare for output
varns={modl,lon,lat};
end
if nargin>3 && flag~=1
defval('toplot',0)
if toplot==1
% Do the plotting in Mollweide
clf ; [d,ch,ph]=plotplm(modl,[],[],1,1);
% Adjust the colorbar to a symmetric one
caxis(round(halverange(modl,75)))
% Make a color scale that is "tomographic"
kelicol
% Put a color bar on
cb=colorbar('hor'); shrink(cb,1.5,1.5);
longticks(cb,2); axes(cb)
xl=xlabel('shear velocity variation from 1-D');
movev(cb,-.05); set(cb,'xaxisl','b')
% Cosmetics
set(cat(1,ch{:}),'linew',1.5);
set(ph,'color','g','linew',1.5)
fig2print(gcf,'portrait')
% Prepare for output
varns={modl,lon,lat,ch,ph,cb,xl};
end
end
% Collect output
varargout=varns(1:nargout);
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.