function [alphaaz,e1,e2,se1,se2]=plot_tensor(lat,lon,e_phiphi,e_thetatheta,...
    e_thetaphi,s_e_phiphi,s_e_thetatheta,s_e_thetaphi,s,opt,gmtflg)
%  function [alphaaz,e1,e2,se1,se2]=plot_tensor(lat,lon,e_phiphi,e_thetatheta,...
%          e_thetaphi,s_e_phiphi,s_e_thetatheta,s_e_thetaphi,s,opt,gmtflg);
%
% plots the sphereical strain tensor in a geographic map view context.
%
% s is scale factor to apply to plot
% opt is option to plot only the tensor elements that are significant
%   (significance now defined as abs(strain)>2sigma)
%
% Note that the tensor elements are plotted with north going up, east to the right.
% alphaaz is the azimuth of the maximum extension direction clockwise from north.
%
% blue is extension; red is contraction
%
% choose option to plot all or just significant tensor elements.
% 1 means to plot all.

if gmtflg   
    fid1=fopen('plot_tensor_e1.txt','w');
    fid2=fopen('plot_tensor_e2.txt','w');
end

if isempty(opt)
    opt=0;
end;

truescale = 0;

A = [e_phiphi -e_thetaphi;
     -e_thetaphi e_thetatheta];
[V,D]=eig(A);
[e1,imax] = max(diag(D));
[e2,imin] = min(diag(D));
vmax = V(:,imax);

alpha = atan2(vmax(2),vmax(1))*180/pi;
alphaaz = 90-alpha; % make it into a geographic azimuth.

% compute uncertainties in the principal strain rates.
B = .5*(e_phiphi - e_thetatheta);
A = .5*(e_thetaphi.^2 + B.^2).^(-.5);
se1= sqrt( ... 
     (s_e_phiphi.*(.5 + A.*B)).^2 + ...
     (s_e_thetatheta.*(.5 - A.*B)).^2 + ...
     (s_e_thetaphi.*2.*A.*e_thetaphi).^2);
se2= sqrt( ... 
     (s_e_phiphi.*(.5 - A.*B)).^2 + ...
     (s_e_thetatheta.*(.5 + A.*B)).^2 + ...
     (s_e_thetaphi.*2.*A.*e_thetaphi).^2);

axpbar = get(gca,'dataaspectratio');
alpha2=alpha*pi/180;

if (abs(e1) > 2*se1 || opt)
    if truescale == 1
        dx = cos(alpha2)*abs(e1*s);
        dy = sin(alpha2)*abs(e1*s);
    else
        dx = cos(alpha2)*abs(e1*s)*axpbar(1);
        dy = sin(alpha2)*abs(e1*s)*axpbar(2);
    end
    if (e1>0)
        sym='b-';
    else
        sym='r-';
    end;
    plot([lon-dx lon+dx],[lat-dy lat+dy],sym,'linewidth',3);

    if gmtflg
        fprintf(fid1,'%9.5f %9.5f\n',[lon-dx lat-dy]);
        fprintf(fid1,'%9.5f %9.5f\n>\n',[lon+dx lat+dy]);
    end;
end

hold on;

if (abs(e2) > 2*se2 || opt)
    if truescale == 1
		dx = cos(alpha2+pi/2)*abs(e2*s);
		dy = sin(alpha2+pi/2)*abs(e2*s);
    else     
        dx = cos(alpha2+pi/2)*abs(e2*s)*axpbar(1);
     	dy = sin(alpha2+pi/2)*abs(e2*s)*axpbar(2);
    end
	if (e2>0)
       sym='b-';
	else
       sym='r-';
	end;
	plot([lon-dx lon+dx],[lat-dy lat+dy],sym,'linewidth',3);
    
    if gmtflg
        fprintf(fid2,'%9.5f %9.5f\n',[lon-dx lat-dy]);
        fprintf(fid2,'%9.5f %9.5f\n>\n',[lon+dx lat+dy]);
    end;
end


if gmtflg
    fclose(fid1);
    fclose(fid2);
end