%Crash Course in using satelliteTrackPlot
%Input arguments:
%   data - four column array
%       Column 1 should be date/time in 'fractional days' (to convert from
%       seconds of the day to fractional days use:
%       fracdays = UTsec./86400 + daynum; 
%       where daynum is the (scalar) day number and UTsec is the
%       (vector/array) of seconds of the day.
%       Column 2 should be magnetic latitude
%       Column 3 should be magnetic local time
%       Column 4 should be the Vy (or whatever other variable you are
%       interested in plotting)
%       *BEST PRACTICE IS TO ONLY HAVE data CONTAIN ONLY ONE PASS AND ONE HEMISPHERE*
%       *PUT THE FUNCTION IN A LOOP IF YOU NEED MANY PASSES*
%   cvar_range - the range for Vy in the format:
%       cvar_range = [minval,maxval];
%   max_height - the height of a line that represents a value of Vy equal to
%       maxval in cvar_range. This is basically a way to control the size
%       of the colored lines that represent values of Vy
%   varargin - is short for variable input argument
%       this is an OPTIONAL input and it is a way to get the plot to behave
%       differently. Some possible options that may be useful are:
%
%           satelliteTrackPlot(data,cvar_range, max_height,'Axes',ax);
%               -makes the function plot to an already existing set of axes
%               with handle ax.
%           
%           satelliteTrackPlot(data,cvar_range, max_height,'Update');
%               -makes the function plot on top of the current axes instead
%               of opening a new figure and drawing a new set of axes on it
%
%           satelliteTrackPlot(data,cvar_range, max_height,'NoTrack');
%               -plots only the colored lines and not the black satellite
%               track
%
%           satelliteTrackPlot(data,cvar_range, max_height,'NoTimeLabel');
%               -Does not plot the time of the start of the data
%
%           satelliteTrackPlot(data,cvar_range, max_height,'MinimumDisplayThreshold',minpercent);
%               -Exludes the lowest minpercent of data, for example if
%               minpercent is .05, the lowest 5% of data is not plotted
%
%           satelliteTrackPlot(data,cvar_range, max_height,'MaximumDisplayThreshold',maxpercent);
%               -Exludes all data greater than maxpercent*maxval from cvar_range, for example if
%               maxpercent is .8 and maxval is 1000, any data with Vy above
%               800 is not plotted 
%
%   Output arguments - the function dumps arrays of line and point objects
%   into the output arguments. Usually you won't need these at all unless
%   you are interested in manually clearing the track and colored lines
%   without closing the figure. To do this call:
%       delete('lines','track');
            
%function [out,track,lines] = satelliteTrackPlot(data, cvar_range, max_height, PlotTitle, varargin)
function [out,track,lines] = satelliteTrackPlot(data, cvar_range, max_height, PlotTitle, initialize,cbarmax,varargin)
            %satelliteTrackPlot generates a plot of satellite data as colored lines
            %along a plot of the satellite track
            %  Data should have 4 columns: time, maglat, mlt, and pf respectively
            % cvar range should be the range of the variable to be plotted
            % using the colorbar. Format should be [rangestart,rangeend]
            % maxheight is the height in radial units of a value of cvar
            % at the top of cvar_range 
            
%  8/2011 Barbara Emery:  Pass in initialize, and assume if 0 is like 'Update' since code does not seem to work with varargin
%           inputaxes = [];
%           initialize = 1;
            trackson = 1;
            returncolorbar = 0;
            trackcolor = [0,0,0];
            singlecolor = [];
            maxvalsonly = 0;
            threshold = 5; %threshold for maxvalsonly 5 means top 5% of largest values retained
            simplemax = 0;
            minthresh = .05;
            maxthresh = 100;
            markersize = 4;
            mark = 'o';
            edgecol = 'none';
            single_point_color = [];
            logfile = [];
            cnormfilter = [];
            minlat = [];
            timelabelson = 1;
            tracksonly=0;
            logfileyear = [];
            track = [];
            lines = [];
            out = [];
            
% NOTE:  nargin always 0 so comment out!           
%           if nargin > 4 
                k=1;
                while k < length(varargin);                    
                    switch varargin{k}
                        case 'Axes'
                            inputaxes = varargin{k+1};
                            k=k+2;
                        case 'ReturnColorBar'
                            returncolorbar = 1;
                            k=k+1;
                        case 'Update'
                            initialize = 0;
                            inputaxes = varargin{k+1};
                            k=k+2;
                        case 'TrackColor'
                            trackcolor = varargin{k+1};
                            k=k+2;
                        case 'Greyscale'
                            singlecolor = 1;
                            k=k+1;
                        case 'MaxValsOnly'
                            maxvalsonly = 1;
                            k=k+1;
                        case 'NoTrack'
                            trackson = 0;
                            k=k+1;
                        case 'NoTimeLabel'
                            timelabelson = 0;
                            k=k+1;
                        case 'SimpleMax'
                            simplemax = 1;
                            k=k+1;
                        case 'Threshold'
                            threshold = varargin{k+1};
                            k=k+2;
                        case 'MinimumDisplayThreshold'
                            minthresh = varargin{k+1};
                            k=k+2;
                        case 'MaximumDisplayThreshold'
                            maxthresh = varargin{k+1};
                            k=k+2;
                        case 'Marker'
                            mark = varargin{k+1};
                            k=k+2;
                        case 'MarkerSize'
                            markersize = varargin{k+1};
                            k=k+2;
                        case 'SingleColor'
                            single_point_color = varargin{k+1};
                            k=k+2;
                        case 'LogFile'
                            logfile = varargin{k+1};
                            k=k+2;
                        case 'LogFileYear'
                            logfileyear = varargin{k+1};
                            k=k+2;
                        case 'NegCnormOnly'
                            cnormfilter = 'cnorm < 0';
                            k=k+1;
                        case 'PosCnormOnly'
                            cnormfilter = 'cnorm > 0';
                            k=k+1;
                        case 'MinLat'
                            minlat = varargin{k+1};
                            k=k+2;
                        case 'TrackOnly'
                            tracksonly = 1;
                            trackson = 1;
                            timelabelson = 0;
                            k=k+1;
                            
                        otherwise
                        error('invalid optional input argument %s',varargin{k});
                    end
                end
%            end

             if(initialize==1)
%            if isempty(inputaxes)
% ScreenSize = [left,bottom,width,height]            
                scsz = get(0,'ScreenSize');                
%                h = figure('Color',[.95,.95,.95],'Position',[.5*(scsz(3)-scsz(4)),0,scsz(4),scsz(4)]);
 h = figure('Color',[.95,.95,.95],'Position',[.5*(scsz(3)-scsz(4)),0,0.5*scsz(4),0.5*scsz(4)]);

%	        h = figure('Color',[.95,.95,.95],'Position',[0.07,0.07,0.90,0.90]);
                ax = axes('Parent',h,'Visible','off');
            
%            else
%                ax = inputaxes;
%           end
%  8/2011 Barbara Emery:  Move title to first only
%            if(initialize)
              draw_polar(ax)
              title(PlotTitle,'fontsize',16)
            end

            hold on
%NS		% This is undefined and will make things quit
		
            %{
            %parse out missing values
            g =  ~isnan(data(:,1)) & ~isnan(data(:,2)) & ~isnan(data(:,3)) & ~isnan(data(:,4));

            data = data(g,:);

            %parse out N polar parts of passes

            gn = data(:,2) >= 45;
            gs = data(:,2) <= -45;
% 9/11 Emery:  Where is NS passed in?
            if(strcmp(NS,'N'))   
                hemidata = data(gn,:); 
            elseif(strcmp(NS,'S'))
                hemidata = data(gs, :);
            end
%NS		% NOT PRINTED OUT???
length(hemidata(:,1))		% NOT PRINTED OUT??
            mlt2rad = @(mlt) (mlt./24).*(2*pi) - pi/2; 

            maglat2dist = @(mlat) -1*(abs(mlat)-90);


            %convert mlat/mlt to distance / angle in radians
            n = 1; passstart(1) = 1;
            for k = 1:length(hemidata(:,1))-1
                deltat = hemidata(k+1,1)-hemidata(k,1);
                if deltat >= 1/24;
                    passend(n) = k;
                    n=n+1;
                    passstart(n) = k+1;
                end
            end

            passend(n) = length(hemidata(:,1));
            %}
            
            %set up conversion functions 
            mlt2rad = @(mlt) (mlt./24).*(2*pi) - pi/2; %angle in radians starting from -y axis, increasing ccw
            rad2mlt = @(rad) (rad + pi/2)./(2*pi).*24; 

            maglat2dist = @(mlat) -1*(abs(mlat)-90); %radial distance
            dist2maglat = @(dist) (90-abs(dist));
   %{
            %set up the colormap
            reds = [linspace(.5,0,30),linspace(0,.8,35),linspace(.8,1,35)];
            blues = [linspace(.5,.3,20),linspace(.3,.8,40),linspace(.8,0,20),zeros(1,20)];
            greens = [linspace(.5,0,20),zeros(1,60),linspace(0,.5,20)];
            
            %{
            reds = [linspace(0,.8,50),linspace(.8,1,50)];
            blues = [linspace(1,.8,50),linspace(.8,0,50)];
            greens = zeros(1,100);
            %}
       
            cmap = [transpose(reds), transpose(greens), transpose(blues)]; %blue to red colormap
            %}
        
            map = colormap(jet(100));
            

            %set the colorgraph properties         
            
            if initialize
                colorbar('Location','SouthOutside','tag','Colorbar');
                cbar = findobj(gcf,'tag','Colorbar');
                
                nticks = 6;
                ticks = linspace(cvar_range(1), cvar_range(2), nticks);
                
                for k = 1:nticks
                    tickstr{k} = num2str(ticks(k),'%.3f');
                end
                
                ticks = (( ticks-cvar_range(1) )./( cvar_range(2)-cvar_range(1) ) )*100+1;
                
                
                set(cbar,'XTick',ticks);
                set(cbar,'XTickLabel',tickstr);
            end
            
            if ~isempty(minlat)
                inbounds = abs(data(:,2)) > minlat;
                data = data(inbounds,:);
            end
            
            hemi = sign(data(:,2));
            data(:,2) = maglat2dist(data(:,2));
            data(:,3) = mlt2rad(data(:,3)); 
            
            if (maxvalsonly && ~tracksonly)
                
               [pos_x, pos_y] = pol2cart(data(:,3),data(:,2));
                
               hold on
               
               if(trackson)
                track = plot(pos_x, pos_y, 'k.','MarkerSize',1,'Color',trackcolor); %plot with the track
               end
               
               if(simplemax)
                   %fprintf('Note that passing MaxValsOnly causes sattrackPlot to retain the largest 5 percent of abs(data) \n');
                   [maxval, maxind] = max(abs(data(:,4)));
                   threshval = (1-threshold/100)*maxval;
                   %fprintf('Threshold: %d, Threshval: %d \n',threshold,threshval);
                   
                   retain = abs(data(:,4)) >= (1-threshold/100)*maxval;
                   data = data(retain,:);
                   
               else
                   
                   [maxval, maxind] = max(abs(data(:,4)));
                   threshval = (1-threshold/100)*maxval;
                   try
                       [maxvals, maxinds] = findpeaks(abs(data(:,4)),'minpeakheight',threshval,'minpeakdistance',10);
                       data = data(maxinds,:);
                   catch ME
                       %fprintf('Findpeaks failed, using simple max algorithm \n');
                       %fprintf('Threshold: %d, Threshval: %d \n',threshold,threshval);
                       
                       retain = abs(data(:,4)) >= (1-threshold/100)*maxval;
                       data = data(retain,:);
                   end
               end

                [pos_x, pos_y] = pol2cart(data(:,3),data(:,2)); %again with revised data set
               
            else

                
                [pos_x, pos_y] = pol2cart(data(:,3),data(:,2));
                
               if(trackson)
               track = plot(pos_x, pos_y, 'k.','MarkerSize',1,'Color',trackcolor);
               end
            
            end

            if ~isempty(logfile)
               fid = fopen(logfile,'a');
               if fid == -1
                   error('Failed to open file for logging: %s',logfile);
               end
            end
            
            if ~isempty(cnormfilter) %check for addition condition set by input args
                    cnormfilter = strcat('~isnan(cnorm) && ',cnormfilter);
                else 
                    cnormfilter = '~isnan(cnorm)';
            end            
            
            if ~tracksonly
                %create an array with cartesian coordinates and values for the variable
                %that will be plotted with colored lines (usually poynting flux)
                plotarr = [pos_x,pos_y,data(:,4)];
                labeled = 0; c = 1; b=1;
                counter = [];

                if(~maxvalsonly)
                    for l = 1:length(plotarr(:,1))
                        val = plotarr(l,3);
                        xval = plotarr(l,1);
                        yval = plotarr(l,2);

                        %normalize val to the colorbar scale or to cbarmax (B. Emery)
                        scale_ex = 0;
		        if (cbarmax>0.)
		          cnorm = abs(val)./cbarmax;
			else
                          cnorm =  (val-cvar_range(1))./(cvar_range(2) - cvar_range(1));
			end

                        if eval(cnormfilter)
                            
                            %{
                            if(~isempty(counter))
                                for c = 1:length(counter)
                                    fprintf('\b');
                                end
                            end

                            counter = sprintf('SatTrackPlot: %d/%d',l,length(plotarr(:,1)));
                            fprintf('SatTrackPlot: %d/%d',l,length(plotarr(:,1)));
                            %}
                            
                            ln = cnorm*max_height*sign(val);

                                colr = ceil(cnorm*100);

                                %account for values that fall outside of the cvar_range
                                if(colr > 100)
                                    colr = 100;
                                elseif(colr < 1);
                                    colr = 1;
                                end

                            if (isempty(singlecolor))

                                color_array = map(colr,:);
                            else

                                color_array = [colr*.9/100+.1,colr*.9/100+.1,colr*.9/100+.1];
                            end

                            ydat = [yval,yval+ln];
                            xdat = [xval,xval];

                            if(~labeled && trackson && timelabelson)
if data(1,3)>12.
                                lab = text(xval,yval,sprintf('%.2f',data(1,1)),'BackgroundColor',[.5,.5,.5],'Color',[1,1,1],'FontSize',8,...
                                'HorizontalAlignment','right','VerticalAlignment','top');
else
 lab=text(xval,yval,sprintf('%.2f',data(1,1)),'BackgroundColor',[.5,.5,.5],'Color',[1,1,1],'FontSize',8,...
                                'HorizontalAlignment','left','VerticalAlignment','top');
end
                                labeled = 1;
                                track = [track,lab];
                            end
                                
                                lines(l) = plot(xdat,ydat,'LineStyle','-','Color',color_array,'LineWidth',.5);
                        end

                    end
                else 
                    for l = 1:length(plotarr(:,1))

                        dat = [data(l,:),hemi(l)];

                        val = plotarr(l,3);
                        xval = plotarr(l,1);
                        yval = plotarr(l,2);

                        %normalize val to the colorbar scale
                        scale_ex = 0;
		        if (cbarmax>0.)
		          cnorm = abs(val)./cbarmax;
			else
                          cnorm =  (val-cvar_range(1))./(cvar_range(2) - cvar_range(1));
			end
                       

                        if(cnorm < minthresh)
                            continue
                        end

                        if abs(val - cvar_range(1)) > abs(cvar_range(2) - cvar_range(1))
                            %fprintf('Value exceeds color scale! Value: %.2f, Cnorm: %.2f \n',val, cnorm);
                            scale_ex = 1;
                        end

                        if eval(cnormfilter)

                            ln = cnorm*max_height*sign(val);

                            %fprintf('PF value: %f, Normalized Value: %f \n',val,cnorm);

                            colr = ceil(cnorm*100);

                            %account for values that fall outside of the cvar_range
                            if(colr > 100)
                                colr = 100;
                            elseif(colr < 1);
                                colr = 1;
                            end

                            if (isempty(singlecolor))

                                color_array = map(colr,:);
                            else

                                color_array = [colr*.9/100+.1,colr*.9/100+.1,colr*.9/100+.1];
                            end

                            ydat = [yval,yval+ln];
                            xdat = [xval,xval];
                            

                            if(~labeled && trackson && timelabelson)
if data(1,3)<12.
                                text(xval,yval,sprintf('%.2f',data(1,1)),'BackgroundColor',[.5,.5,.5],'Color',[1,1,1],'FontSize',8,...
                                'HorizontalAlignment','right','VerticalAlignment','top');
else
 text(xval,yval,sprintf('%.2f',data(1,1)),'BackgroundColor',[.5,.5,.5],'Color',[1,1,1],'FontSize',8,...
                                'HorizontalAlignment','left','VerticalAlignment','top');
end
                                labeled = 1;
                                track = [track,lab];
                            end

                            if abs(cnorm) > minthresh && abs(cnorm) < maxthresh
                                
                                if scale_ex 
                                    fprintf('Large value plotted, cnorm was %.2, color was R: %.2f G: %2f B%.2f \n',cnorm,color_array(1),color_array(2),color_array(3));
                                end

                                if( isempty(single_point_color) )
                                %fprintf('Point Value: %.3f, Cnorm: %.3f min:%.2f max:%2f \n',val,cnorm,minthresh,maxthresh);
                                    lines(l) = plot(xdat(1),ydat(1),'Marker',mark,'Color',color_array,...
                                        'MarkerFaceColor',color_array,'MarkerSize',markersize,'LineWidth',3);

                                    if (~isempty(logfile))
                                        if isempty(logfileyear)
                                            fprintf(fid,'Time: %f, Mlat: %f, MLT: %f, Sz: %f, Sznorm: %f\n',...
                                            dat(1),dat(5)*dist2maglat(dat(2)),rad2mlt(dat(3)),dat(4),cnorm);
                                        else 
                                            fprintf(fid,'%s, Time: %f, Mlat: %f, MLT: %f, Sz: %f, Sznorm: %f\n',...
                                            logfileyear,dat(1),dat(5)*dist2maglat(dat(2)),rad2mlt(dat(3)),dat(4),cnorm);
                                        end
                                    end
                                else
                                     lines(l) = plot(xdat(1),ydat(1),'Marker',mark,'Color',single_point_color,...
                                    'MarkerFaceColor',single_point_color,'MarkerSize',markersize,'LineWidth',3);

                                    if (~isempty(logfile))
                                        if isempty(logfileyear)
                                            fprintf(fid,'Time: %f, Mlat: %f, MLT: %f, Sz: %f, Sznorm: %f\n',...
                                            dat(1),dat(5)*dist2maglat(dat(2)),rad2mlt(dat(3)),dat(4),cnorm);
                                        else 
                                            fprintf(fid,'%s: Time: %f, Mlat: %f, MLT: %f, Sz: %f, Sznorm: %f\n',...
                                            logfileyear,dat(1),dat(5)*dist2maglat(dat(2)),rad2mlt(dat(3)),dat(4),cnorm);
                                        end
                                    end
                                end
                                
                            end
                    end 
                end
                end
                
            end

      
            
           
            %text(xval,yval,sprintf('%0.2f',data(end,1)-floor(data(end,1))),'BackgroundColor',[.5,.5,.5],'Color',[1,1,1],'FontSize',8,...
                   % 'HorizontalAlignment','right','VerticalAlignment','bottom');
            
            
            if returncolorbar && initialize
                out = cbar;
            end
            
            if ~isempty(logfile)
                fclose(fid);
            end
        end  

function draw_polar(ax)
            %This code written by Adela Anghel; copied from Data_Plotter program
            axes(ax)
            set(gca,'Color',[.95,.95,.95]);
            FS = 15;
            FS1 = 15;
            color_grid = [0.10 0.10 0.10];
            phi = linspace(0,2*pi,10000); 
            linestyle = ':';
            minlat = 50;

            switch minlat
                case 60
                    
                    plot(30*cos(phi), 30*sin(phi),'LineStyle','-','Color',color_grid,...
                         'LineWidth',1);
                    axis equal; hold on; axis off
                    plot(20*cos(phi), 20*sin(phi),linestyle,'Color',color_grid);
                    plot(10*cos(phi), 10*sin(phi),linestyle,'Color',color_grid);
                
                case 50
                    
                    plot(40*cos(phi), 40*sin(phi),'LineStyle','-','Color',color_grid,...
                         'LineWidth',1);
                    axis equal; hold on; axis off
                    plot(30*cos(phi), 30*sin(phi),linestyle,'Color',color_grid);
                    plot(20*cos(phi), 20*sin(phi),linestyle,'Color',color_grid);
                    plot(10*cos(phi), 10*sin(phi),linestyle,'Color',color_grid);

                case 40
            
                    plot(50*cos(phi), 50*sin(phi),'LineStyle','-','Color',color_grid,...
                         'LineWidth',1);
                    axis equal; hold on; axis off
                    plot(40*cos(phi), 40*sin(phi),linestyle,'Color',color_grid); 
                    plot(30*cos(phi), 30*sin(phi),linestyle,'Color',color_grid);
                    plot(20*cos(phi), 20*sin(phi),linestyle,'Color',color_grid);
                    plot(10*cos(phi), 10*sin(phi),linestyle,'Color',color_grid);
                    
                otherwise
                    
                    error('Minlat is not one of the supported values: 60, 50, 40');
                    
             end
                    r_text = 90-minlat+2; th_text = 3*pi/2;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-1,y,'0','FontSize',FS,'Color','r')

                    r_text = 90-minlat+2; th_text = 7*pi/4;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-1,y,'3','FontSize',FS,'Color','r')

                    r_text = 90-minlat+2; th_text = 0*pi/2;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-1,y,'6','FontSize',FS,'Color','r')

                    r_text = 90-minlat+2; th_text = pi/4;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-1,y,'9','FontSize',FS,'Color','r')

                    r_text = 90-minlat+2; th_text = 1*pi/2;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-2,y,'12','FontSize',FS,'Color','r')

                    r_text = 90-minlat+2; th_text = 3*pi/4;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-2,y,'15','FontSize',FS,'Color','r')

                    r_text = 90-minlat+2; th_text = 2*pi/2;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-1,y,'18','FontSize',FS,'Color','r')

                    r_text = 90-minlat+2; th_text = 5*pi/4;
                    [x,y] = pol2cart(th_text,r_text);
                    text(x-2,y,'21','FontSize',FS,'Color','r')
            
            % line([0 0],[-50 50],'Color',[0.7 0.7 0.7],'LineWidth',1.5); 
            % line([-50 50],[0 0],'Color',[0.7 0.7 0.7],'LineWidth',1.5)
            for i = 1:24
                th = (i-1)*pi/12;
                r_min = 10; 
                r_max = 90-minlat;
                
                line([r_min*cos(th) r_max*cos(th)],[r_min*sin(th) r_max*sin(th)],...
                     'Color',color_grid,'LineStyle',linestyle,'LineWidth',1)
            end

            text( 4,-5,'80^o','FontSize',FS1,'Color','r');
            text(14,-5,'70^o','FontSize',FS1,'Color','r');
            text(24,-5,'60^o','FontSize',FS1,'Color','r');
            
            if(minlat < 60)
            text(34,-5,'50^o','FontSize',FS1,'Color','r');
            end
            
            if(minlat < 50)
            text(44,-5,'40^o','FontSize',FS1,'Color','r');
            end
      
            xlim([-(r_max+5),(r_max+5)]); ylim([-(r_max+5),(r_max+5)]);
        end
