function [varargout]=plotraypaths(paths,shells,cmap,ax)
%PLOTRAYPATHS    Plot raypaths in 3D
%
%    Usage:  plotraypaths(paths)
%            plotraypaths(paths,shells)
%            plotraypaths(paths,shells,ax)
%            ax=plotraypaths(...)
%
%    Description:
%     PLOTRAYPATHS(PATHS) draws the raypaths in the TauP struct PATHS
%     (generated by TAUPPATH or GETRAYPATHS) in a 3D Earth plot.  The
%     surface is covered with an topography image and the transparency is
%     adjusted so that the oceans are highly transparent.  Shells
%     corresponding to the 410, 670, CMB, & ICB are also plotted to help
%     visualize where the rays pierce those depths.
%
%     PLOTRAYPATHS(PATHS,SHELLS) sets the depth shells plotted within the
%     Earth.  SHELLS is a Nx3 cell array of {depth1 color1 alpha1; ...}
%     where depth is in km, color may be a rgb triplet or colorname and
%     alpha is the transparency of the shell (0 to 1).  N is the number of
%     shells.  The default is: { 400   'b' .2;
%                                670   'g' .2;
%                               2891   'o' .7;
%                               5149.5 'r' .7}.
%
%     PLOTRAYPATHS(PATHS,SHELLS,CMAP) alters the raypath colormap.  The
%     default is the 'hsv'.  CMAP may be a string/handle to a colormap
%     function or a NRPx3 array where NRP is the number of raypaths in
%     PATHS.
%
%     PLOTRAYPATHS(PATHS,SHELLS,CMAP,AX) sets which axis to plot the
%     raypaths in.  This is useful for subplots and guis.
%
%     AX=PLOTRAYPATHS(...) outputs the axis used for the plot.
%
%    Notes:
%
%    Examples:
%     % Plot some seismic phases between a random station & event:
%     plotraypaths(tauppath('s',randlatlon,'e',randlatlon,'z',700*rand))
%
%     % Get several phases for a SEIZMO array dataset and plot them up:
%     Ppaths=getraypaths(data,'P','prem');
%     PcPpaths=getraypaths(data,'PcP','prem');
%     PPpaths=getraypaths(data,'PP','prem');
%     plotraypaths([Ppaths PcPpaths PPpaths]);
%
%    See also: GETRAYPATHS, TAUPPATH, PLOT_TAUPPATH

%     Version History:
%        June  1, 2010 - initial version
%        June  3, 2010 - handle empty path
%        Jan.  4, 2011 - added example, improved plot calls
%        Feb. 27, 2012 - update for tauppath & getraypaths changes
%
%     Written by Garrett Euler (ggeuler at wustl dot edu)
%     Last Updated Feb. 27, 2012 at 02:15 GMT

% todo:

% check nargin
error(nargchk(1,4,nargin));

% check raypath struct
test=tauppath('ph','P','ev',[0 0],'st',[0 10]);
if(~isstruct(paths) || any(~ismember(fieldnames(paths),fieldnames(test))))
    error('seizmo:plotraypaths:badStruct',...
        'PATHS does not appear to be a valid raypath struct!');
elseif(any(~ismember(fieldnames(paths(1).path),fieldnames(test(1).path))))
    error('seizmo:plotraypaths:badStruct',...
        'PATHS does not appear to be a valid raypath struct!');
elseif(any(~ismember({'latitude' 'longitude'},fieldnames(paths(1).path))))
    error('seizmo:plotraypaths:badStruct',...
        'Latitude & Longitude are required path fields!');
end

% number of raypaths
nrp=numel(paths);

% default/check shells
if(nargin==1 || isempty(shells))
    shells={400 'b' .2; 670 'g' .2; 2891 'o' .7; 5149.5 'r' .7};
end
if(~iscell(shells) || size(shells,2)~=3)
    error('seizmo:plotraypaths:badSHELL',...
        'SHELLS must be a Nx3 cell array!');
elseif(~all(cellfun('isreal',shells(:,1))) ...
        || ~all(cellfun('prodofsize',shells(:,1))==1))
    error('seizmo:plotraypaths:badSHELL',...
        'Depth values in SHELLS must be real-valued scalars!');
elseif(~all(cellfun('isreal',shells(:,3))) ...
        || ~all(cellfun('prodofsize',shells(:,3))==1) ...
        || any([shells{:,3}]<0 | [shells{:,3}]>1))
    error('seizmo:plotraypaths:badSHELL',...
        'Alpha values in SHELLS must be scalars within 0 & 1!');
end

% check shell colors
nsh=size(shells,1);
isc=cellfun('isclass',shells(:,2),'char');
shells(isc,2)=mat2cell(name2rgb(shells(isc,2)),ones(sum(isc),1));
if(~all(cellfun('prodofsize',shells(:,2))==3))
    error('seizmo:plotraypaths:badSHELL',...
        'Color values in SHELLS must be rgb triplets or colornames!');
end

% check colormap
if(nargin<3 || isempty(cmap)); cmap=hsv(nrp); end
if(isreal(cmap) && (~isequal(size(cmap),[nrp 3]) ...
        || any(cmap(:)<0 | cmap(:)>1)))
    error('seizmo:plotraypaths:badCMAP',...
        'CMAP must be a Nx3 rgb array or a str/handle of a colormap!');
elseif(isreal(cmap))
    % good to go
elseif(ischar(cmap))
    cmap=str2func(cmap);
    cmap=cmap(nrp);
elseif(isa(cmap,'function_handle'))
    cmap=cmap(nrp);
else
    error('seizmo:plotraypaths:badCMAP',...
        'CMAP must be a Nx3 rgb array or a str/handle of a colormap!');
end

% check axis
if(nargin<4 || isempty(ax) || ~isscalar(ax) || ~isreal(ax) ...
        || ~ishandle(ax) || ~strcmp('axes',get(ax,'type')))
    fh=figure('color','k','defaulttextcolor','w',...
        'defaultaxesxcolor','w','defaultaxesycolor','w',...
        'name','3D Ray Paths','tag','plotraypaths');
    ax=axes('parent',fh);
else
    fh=get(ax,'parent');
    
    % clear any existing objects in the axes
    h=get(ax,'children'); delete(h);
    h=findobj(get(fh,'children'),'peer',ax); delete(h);
end

% constants
Re=6371;

% loop over each raypath
rph=nan(1,nrp);
cnt=0;
for i=1:nrp
    % skip if no points
    if(isempty(paths(i).path.depth)); continue; end
    cnt=cnt+1;
    
    % convert lat/lon/depth to x/y/z
    [x,y,z]=geocentric2xyz(paths(i).path.latitude,...
        paths(i).path.longitude,paths(i).path.depth,Re);
    
    % plot path
    rph(i)=plot3(-x,-y,z,'color',cmap(i,:),'parent',ax);
    if(cnt==1); hold(ax,'on'); end
    
    % plot start/end points
    [evx,evy,evz]=geocentric2xyz(paths(i).event(1),...
        paths(i).event(2),0,Re);
    [stx,sty,stz]=geocentric2xyz(paths(i).station(1),...
        paths(i).station(2),0,Re);
    if(cnt==1)
        h1=plot3(-evx,-evy,evz,'yp','parent',ax,...
            'markersize',10,'markerfacecolor','y');
        h2=plot3(-stx,-sty,stz,'r^','parent',ax,...
            'markersize',4,'markerfacecolor','r');
    else
        plot3(-evx,-evy,evz,'yp','parent',ax,...
            'markersize',10,'markerfacecolor','y');
        plot3(-stx,-sty,stz,'r^','parent',ax,...
            'markersize',4,'markerfacecolor','r');
    end
end

% plot x y z axes
%[x,y,z]=geocentric2xyz(0,0,[-629 6371],6371);
%plot3(-x,-y,z,'g','linewidth',2,'parent',ax);
%[x,y,z]=geocentric2xyz(0,90,[-629 6371],6371);
%plot3(-x,-y,z,'b','linewidth',2,'parent',ax);
%[x,y,z]=geocentric2xyz(90,0,[-629 6371],6371);
%plot3(-x,-y,z,'r','linewidth',2,'parent',ax);

% plot cmb/icb
[x,y,z]=sphere(45);
for i=1:nsh
    r=Re-shells{i,1};
    surface(r*x,r*y,r*z,...
        'parent',ax,...
        'facecolor',shells{i,2},...
        'facealpha',shells{i,3},...
        'backfacelighting','unlit',...
        'edgecolor','none');
end

% load topography
tmp=load('topo.mat');
if(~isfield(tmp,'topo') || ~isreal(tmp.topo) || ndims(tmp.topo)~=2 ...
        || size(tmp.topo,1)<8 || size(tmp.topo,2)<16)
    error('seizmo:plotraypaths:badTOPO',...
        'The file topo.mat appears corrupted!');
end

% tranparent globe with continents
s=surface(Re*x,Re*y,Re*z,...
    'facecolor','texturemap','cdata',tmp.topo,'parent',ax);
set(s,'edgecolor','none','facealpha','texture','alphadata',tmp.topo);
set(s,'backfacelighting','unlit');
colormap(ax,sealand);
alpha(s,'direct');
alphamap(fh,[.1;.8])
axis(ax,'off','vis3d');
view(ax,[-90 0])
rotate3d(ax,'on');
zoom(ax,1.5);

% legend
lh=legend(ax,[h1 h2],{'event' 'station'},...
    'location','westoutside','textcolor','w');
set(lh,'color','none','edgecolor','w',...
    'fontsize',6,'interpreter','none');
hold(ax,'off');

% output
if(nargout); varargout{1}=ax; end

end
