% TA_sphere_demo
% Demo code showing the basics of how to create a spherical
% projection of the Earth and how to plot symnbols. This can be extended to
% create a movie, using the same straategy as in the code
% TA_wavefield_demo.m
%
% Bob Woodward
% Incorporated Research Institutions for Seismology / IRIS
% woodward@iris.edu
%
clf
%
% Read a topo database - this sample is extracted from ETOPO1 and saved in
% matlab format.
% The array 'values" will contain a lat-lon grid, with a elev value at
% each array element. The array "topomap1" is a 64x3 colormap array,
% though we are not using it here - we are using the topomap3 colormap
% instead.
%
load('topohires.mat','values','topomap1');
load topomap3
colormap(topomap3);
%
% Read in political boundaries and transform to spherical coordinates.
% You can get coastlines and political boundaries from NOAA at
% http://www.ngdc.noaa.gov/mgg_coastline/
% Use the pull down menu labeled "coastline database" to select either a
% coastline or political boundary database. Note the option for matlab
% format output.
%
political = importdata('political_US.dat');
political = political * pi / 180.0;
[xpol ypol zpol] = sph2cart(political(:,1), political(:,2), 1.0);
%
% create a sphere of unit radius, using a built-in function
[x,y,z] = sphere(180);
%
% set up the plotting
% Most of these values have been discovered through trial and error, or
% are based on values used in matlab-provided demos (in the matlab help
% menu try searching for the "Earth's Topography" demo).
%
cla reset
axis square off
props.AmbientStrength = 0.1;
props.DiffuseStrength = 1;
props.SpecularColorReflectance = .5;
props.SpecularExponent = 20;
props.SpecularStrength = 1;
props.FaceColor        = 'texture';
props.EdgeColor        = 'none';
props.FaceLighting     = 'phong';
caxis([-7000. 7000.]); % limits of the elev/bathymetry data we are plotting
props.Cdata = values;  % identify the values for the sphere's surface
%
% Plot the sphere, with topo / bathymetry
surface(x,y,z,props);
%
% set up the light source position
theta = 205;
phi   = 34;
r     = 2;
[lx ly lz] = sph2cart(theta*pi/180, phi*pi/180, r);
light('position',[lx ly lz]);
%
% define the viewing azimuth and elevation
view(-22,40)
%
% this value effectively sets the zoom level (smaller number zooms in more)
set(gca,'CameraViewAngle',6.0);
%
% set up some data to plot
%
% first, a fictitious earthquake location
eqlat = 36.1;
eqlon = -113.0;
%
% now get coordinates of a small circle and a big circle on surface of the
% sphere
i = 0;
for az = 0 : 30 : 360.0     % step in 30 degreee azimuth increments
    i=i+1;
    [s_circ(i,1) s_circ(i,2)] = givloc(eqlat,eqlon, 1.0,az); % a 1 degree radius circle
    [b_circ(i,1) b_circ(i,2)] = givloc(eqlat,eqlon,10.0,az); % a 10 degree radius circle
end;
%
% Now, experiment with different ways to plot on the sphere.
%
hold on
%
% plot political boundaries on surface of sphere
plot3(xpol,ypol,zpol,'k');
%
% Create a cyan patch by filling the small, 1 degree radius circle.
% Note that in the converstion to cartesion coordinats we use a radius
% slightly greater than 1. Since the filled circle will be a flat surface
% we want it to be slightly above the sphere surface - otherwise it will
% intersect the surface.
s_circ = s_circ * pi / 180.0;
[xx yy zz] = sph2cart(s_circ(:,2), s_circ(:,1), 1.01);
patch(xx,yy,zz,'c');
%
% plot the 10 degree radius circle (polygon) - in yellow
b_circ = b_circ * pi / 180.0;
[xx yy zz] = sph2cart(b_circ(:,2), b_circ(:,1), 1.0);
plot3(xx,yy,zz,'y');
%
% plot a red sphere at 5 degree latitude offset from red sphere
% area of red sphere is 100 square points
[x2,y2,z2] = sph2cart(eqlon*pi/180, (eqlat+5)*pi/180, 1.0);
scatter3(x2,y2,z2,100,[1 0 0 ],'filled');
%
% Try writing some text on sphere - note that we set the radius to a 
% little above the sphere surface - otherwise the text pokes through 
% the surface (e.g., experiment with setting the radius to 1.0).
[x y z] = sph2cart(-128*pi/180, 46*pi/180, 1.1);
text(x,y,z,'sample text','FontSize',20,'Color','white','FontWeight','bold');

