function SPIE2009_4
% SPIE2009_4
%
% A plot of the L=60 expanded gravitational potential perturbations due
% to some large end-member yearthquakes.
%
% Last modified by fjsimons-at-alum.mit.edu, 06/03/2010

disp('Until normalization is fixed we will not have a good plot')

% Produce fictitious earthquake at the equator
defval('deplonlat',[31 95 5])
defval('degres',0.5)
defval('Lup',60)

% Make the moment tensors, see Dahlen & Tromp Table 5.1
M{1}=[ 1  0 -1  0 0 0 ]/sqrt(2); 
M{2}=[ 0  0  0  0 1 0 ]/sqrt(2); 
M{3}=[ 0  1 -1  0 0 0 ]/sqrt(2); 
M{4}=[ 1 -1  0  0 0 0 ]/sqrt(2); 
M{5}=[ 0  0  0  1 0 0 ]/sqrt(2); 
M{6}=[ 0  0  0  0 0 1 ]/-sqrt(2); 

% Produce the first of the potential perturbations, collect eigenfunctions
[p{1},pr{1},Sstrain]=eqpotential(M{1},deplonlat,Lup);

% Produce the rest, reuse eigenfunctions
for ind=2:length(M)
  [p{ind},pr{ind},Sstrain]=eqpotential(M{ind},deplonlat,Lup,Sstrain);
end

% Then do the plotting, for speed only expand locally
TH=10;
c11=[deplonlat(2)-2*TH deplonlat(3)+2*TH];
cmn=[deplonlat(2)+2*TH deplonlat(3)-2*TH];

clf
[ah,ha,H]=krijetem(subnum(3,2));

% Filter ever so slightly before expanding? Don't bother
[r,lon,lat,Plm]=plm2xyz(plmfilt(pr{1},Lup-2,5),degres,[c11 cmn]);
axes(ha(1))
% Plot rotated function on the map
imagefnan(c11,cmn,setnans(r),...
	  kelicol,halverange(r));
[ac,pc(1)]=plotcont(c11+[-2 2],cmn+[2 -2]);
t(1)=title(sprintf('%s = [%i %i %i %i %i %i]/%s','$\mathbf{M}$',...
		   sqrt(2)*M{1},'$\sqrt{2}$'),'interpreter','latex');

for ind=2:length(M)
  r=plm2xyz(plmfilt(pr{ind},Lup-2,5),degres,[c11 cmn],[],[],Plm);
  axes(ha(ind))
  imagefnan(c11,cmn,setnans(r),...
	    kelicol,halverange(r));
  [ac,pc(ind)]=plotcont(c11+[-2 2],cmn+[2 -2]);
  t(ind)=title(sprintf('%s = [%i %i %i %i %i %i]/%s','$\mathbf{M}$',...
		   sqrt(2)*M{ind},'$\sqrt{2}$'),'interpreter','latex');
end

% Cosmetics
nolabels(ah(1:4),1)
nolabels(ha(4:6),2)
set(pc,'linew',0.5)
longticks(ah)
serre(H,1/2,'across')
set(t,'FontS',12)
%moveh(ha(1:3),0.1)
set(ah,'xlim',[c11(1) cmn(1)],'ylim',[cmn(2) c11(2)],...
       'xtick',[80:10:110],'ytick',[-10:10:20])
[bhl,thl]=label(ah,'ll',[],[],[],[],[],[],[]);
deggies(ah)
fig2print(gcf,'tall')
figdisp([],[],[],0)