% =============================================================================
% ============================== PulseGraph ===================================
% =============================================================================
%      ________________________________________________________________
%     |   Approximate the lighting channel by a graph using VHF pulse  |
%     |   locations as vertices in the graph. The pulse data should    |
%     |   be in the format generated by New Mexico Tech's lma_analysis |
%     |   program. The graph is obtained by connecting together        |
%     |   closest neighboring pulse. As an option, the pulse graph and |
%     |   VHF pulses can be plotted. Both stationary 3D plots and      |
%     |   animations can be generated.                                 |
%     |                                                                |
%     |                  Version 1.1  (January 31, 2011)               |
%     |                  Version 1.2  (May 2, 2011)                    |
%     |                  Version 2.0  (June 19, 2012)                  |
%     |                  Version 2.1  (March 5, 2013)                  |
%     |                                                                |
%     |                       William W. Hager                         |
%     |                         hager@ufl.edu                          |
%     |                 http://www.math.ufl.edu/~hager                 |
%     |                   Department of Mathematics                    |
%     |                     University of Florida                      |
%     |                 Gainesville, Florida 32611 USA                 |
%     |                      352-392-0281 x 244                        |
%     |                                                                |
%     |                 Copyright by William W. Hager                  |
%     |                                                                |
%     |This program is free software; you can redistribute it and/or   |
%     |modify it under the terms of the GNU General Public License as  |
%     |published by the Free Software Foundation; either version 2 of  |
%     |the License, or (at your option) any later version.             |
%     |This program is distributed in the hope that it will be useful, |
%     |but WITHOUT ANY WARRANTY; without even the implied warranty of  |
%     |MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the   |
%     |GNU General Public License for more details.                    |
%     |                                                                |
%     |You should have received a copy of the GNU General Public       |
%     |License along with this program; if not, write to the Free      |
%     |Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
%     |MA  02110-1301  USA                                             |
%     |________________________________________________________________|
%

% ========================= INPUT =============================================
% Input:
%     LMAfilename
%     Optional parameters: 'plot', 'rotate', 'noaxes', 'avi', 'parm', 'quiet'

% This code is in the form of a MATLAB function. The first input argument
% is required while the remaining input arguments are optional. The
% optional arguments can be given in any order.

% 1. LMAfilename - (input) In quotes give either the full or relative
%              path name of the file containing the LMA data. The file
%              should contain only the LMA data for the flash, not the
%              header information. A flash file generated by the
%              PartitionStorm code is ideal input. You can also directly
%              input the LMA time, latitude, longitude, and altitude by
%              using 'input' as the LMAfilename. In this case, the
%              next argument is an array with 4 columns containing the
%              time, latitude, longitude, and altitude respectively.

% 2. 'plot'  - (optional input) Generate a 3D stationary plot of
%              the pulses and the pulse graph.

% 3. 'rotate'- (optional input) Generate a 3D plot of the pulses and the
%              pulse graph and rotate the plot by 360 degrees.

% 4. 'noaxes'- (optional input) No axes appear in the plot.

% 5. 'avi'   - (optional input) Generate an avi video, the next argument
%              after the keyword is the name of the avi file to store the video

% 6. 'parm'  - (optional input) Print the parameters.

% 7. 'quiet' - (optional input) Eliminate the informative messages.

% ========================= OUTPUT ============================================
% Output:
%     EdgeList
%     Optional: Vertex1, Vertex2

% 1. EdgeList- (output) An array with two rows. The i-th column of the array
%              contains a pair of integers corresponding to the pulse
%              numbers in the LMA file that are connected together to
%              in forming the pulse graph

% 2. Vertex1 and Vertex2 - (optional output) Arrays with 3 rows,
%              the i-th column of the arrays contains the [x, y, z]
%              coordinates of pulse pairs that are connected together
%              to form the pulse graph. The coordinates are relative
%              to a point of interest (POI) specified in the paramater
%              file PulseGraphParm.m.

% ========================= ALGORITHM =========================================
% Description of the algorithm for generating the pulse graph of a flash:
% Neighboring pulses are connected together to form line segments.
% Line segments are connected together to obtain channel segments.
% Channel segments are joined to obtain larger segments. The
% growth process continues, subject to the thresholds given in
% the file PulseGraphParm.m. We only join nodes that are relatively
% close together, and a node far from the others is discarded.
% =============================================================================

function [EdgeList, Vertex1, Vertex2] = PulseGraph_new (LMAfilename, varargin)
% function [EdgeList, Vertex1, Vertex2, ndangling] = PulseGraph (LMAfilename, varargin)

% Get the path
% fp = fopen ('AnalysisPathName') ;
% if ( fp < 0 )
%     fprintf (1, '\nBefore running PulseGraph, the setup program should be\n') ;
%     fprintf (1, 'used to record the path for the StormAnalysis directory.\n') ;
%     fprintf (1, 'Change to the StormAnalysis directory,\n') ;
%     fprintf (1, 'start MATLAB, and type setup at the prompt.\n\n') ;
%     error(' ') ;
% end
% PathName = fgetl (fp) ;
PathName = '/home/jeff/Documents/MATLAB/StormAnalysis-2.1';
addpath (strcat (PathName, '/Parm')) ;
addpath (strcat (PathName, '/Topo')) ;

if ( strcmp(LMAfilename, 'input') )
    LMA = varargin{1} ;
    LMAt = LMA (:, 1) ;
    LMAlat = LMA (:, 2) ;
    LMAlon = LMA (:, 3) ;
    LMAalt = LMA (:, 4) ;
    i = 2 ;
else
    arg_offset = 1 ;
    flma = fopen (LMAfilename, 'r') ;
    if ( flma < 0 )
        fprintf (1, 'The LMA filename %s that was given does not exist\n', ...
                 LMAfilename) ;
        error ('stop') ;
    end

    % Read the LMA data
    lma_cell = textscan (flma, '%f%f%f%f%f%f%s') ;
    fclose (flma) ;
    LMAt = lma_cell{1} ;
    LMAlat = lma_cell{2} ;
    LMAlon = lma_cell{3} ;
    LMAalt = lma_cell{4} ;
    i = 1 ;
end

PlotGraph = 0 ;
Rotate = 0 ;
NoAxes = 0 ;
Avi = 0 ;
Parm = 0 ;
Quiet = 0 ;
if ( nargin > 1 )
    while ( i <= nargin-1 )
        found = 0 ;
        if ( strcmp('plot', varargin{i}) )
            PlotGraph = 1 ;
            found = 1 ;
        elseif ( strcmp('rotate', varargin{i}) )
            Rotate = 1 ;
            found = 1 ;
        elseif ( strcmp('noaxes', varargin{i}) )
            NoAxes = 1 ;
            found = 1 ;
        elseif ( strcmp('parm', varargin{i}) )
            Parm = 1 ;
            found = 1 ;
        elseif ( strcmp('quiet', varargin{i}) )
            Quiet = 1 ;
            found = 1 ;
        elseif ( strcmp('avi', varargin{i}) )
            Avi = 1 ;
            i = i + 1 ;
            AviFilename = char(varargin{i}) ;
            found = 1 ;
        end
        if ( found == 0 )
            fprintf (1, 'Input parameter "%s" is invalid for ', ...
                         char(varargin{i})) ;
            fprintf (1, 'the PulseGraph code\n') ;
            error ('stop') ;
        end
        i = i + 1 ;
    end
end

% The parameters associated with the algorithm, topographic map,
% point of interest, and plot are contained in the file
% PulseGraphParm.m, which is now read.

PulseGraphParm

% Print the parameters
if ( Parm )
    fprintf (1, '==========================================================') ;
    fprintf (1, '===================================\n') ;
    fprintf (1, ...
         'Parmeters values [time (sec), distance (m), angles (degrees)]:\n\n') ;
    fprintf (1, '\n1. Algorithm parameters.\n') ;
    fprintf (1, '  node_connect ................ %7i\n', node_connect) ;
    fprintf (1, '  channel_connect ............. %7i\n', channel_connect) ;
    fprintf (1, '  RemoveDangle    ............. %7i\n', RemoveDangle) ;
    fprintf (1, '  DanglingChannel ............. %7i\n', DanglingChannel) ;
    fprintf (1, '\n2. Topographic map (%s).\n', ...
                     TopoDescription) ;
    fprintf (1, '  lat_north  .................. %13.5f\n', lat_north) ;
    fprintf (1, '  lat_south  .................. %13.5f\n', lat_south) ;
    fprintf (1, '  lon_east   .................. %13.5f\n', lon_east) ;
    fprintf (1, '  lon_west   .................. %13.5f\n', lon_west) ;
    fprintf (1, '  map_dimx   .................. %7i\n', map_dimx) ;
    fprintf (1, '  map_dimy   .................. %7i\n', map_dimy) ;
    fprintf (1, '\n3. Point of interest.\n') ;
    fprintf (1, '  POI_lat    .................. %13.5f\n', POI_lat) ;
    fprintf (1, '  POI_lon    .................. %13.5f\n', POI_lon) ;
    fprintf (1, '  POI_alt    .................. %13.5f\n', POI_alt) ;
    fprintf (1, '  Flat       .................. %13.5f\n', Flat) ;
    fprintf (1, '  Flon       .................. %13.5f\n', Flon) ;
    fprintf (1, '\n4. Plot parameters.\n') ;
    fprintf (1, '  Camera_dr  .................. %13.5f\n', Camera_dr) ;
    fprintf (1, '  Camera_dz  .................. %13.5f\n', Camera_dz) ;
    fprintf (1, '  Camera_initial_angle ........ %13.5f\n', ...
                                                         Camera_initial_angle) ;
    fprintf (1, '  Camera_frames ............... %13.5f\n', Camera_frames) ;
    fprintf (1, '  Camera_fps    ............... %13.5f\n', Camera_fps) ;
    fprintf (1, '  FigureHandle  ............... %13.5f\n', FigureHandle) ;
    fprintf (1, '==========================================================') ;
    fprintf (1, '===================================\n') ;
end

nlma = length(LMAt) ;
if ( ~Quiet )
   fprintf (1, 'LMA data read\n') ;
   fprintf (1, 'Initial number LMA pulses: %i\n', nlma) ;
end

% Convert from latitude and longitude to x and y relative to
% point of interest (POI)
LMA = zeros (3, nlma) ;
LMA (1, :) = (LMAlon-POI_lon)*Flon ;
LMA (2, :) = (LMAlat-POI_lat)*Flat ;
LMA (3, :) = LMAalt' ;

nisolated = 0 ;
LMA_ones = ones (1, nlma) ;
% Initial pass, find all the isolated pulses and find the average distance
% between nonisolated pulses.
dist = zeros (nlma, 1) ;
avedist = 0 ;
n = 0 ;
for j = 1:nlma
    R = LMA - LMA (:, j)*LMA_ones ;
    Rnorm = sum (R.^2) ;
    Rnorm (j) = inf ;
    [rmin k] = min (Rnorm) ;
    rmin = sqrt(rmin) ;
    if ( rmin > node_connect ) % skip this pulse
        LMA_ones (j) = 0 ;
        nisolated = nisolated + 1 ;
        dist (j) = 0 ;
    else
        avedist = avedist + rmin ;
        dist (j) = rmin ;
        n = n + 1 ;
    end
end
avedist = avedist/n ;
if ( ~Quiet )
    fprintf (1, 'Average distance between nonisolated pulses: %i\n', avedist) ;
end

% Also delete pulses that are relatively far from the other pulses.
cutoff = IsoFactor*avedist ;
for j = 1:nlma
    if ( dist (j) > cutoff )
        LMA_ones (j) = 0 ;
        nisolated = nisolated + 1 ;
    end
end
    
if ( ~Quiet )
    fprintf (1, 'Number of  isolated nodes: %i\n', nisolated) ;
end

% remove all isolated nodes

LMA_keep = find (LMA_ones) ;
LMA = LMA (:, LMA_keep) ;
nlma = size (LMA, 2) ;
fprintf (1, 'Number of LMA pulses (after removing isolates): %i\n', nlma) ;
LMA_ones = ones (1, nlma) ;
EdgeList = zeros (2, nlma) ;
Sets = zeros (nlma, 1) ;
Label = -ones (nlma, 1) ; % initialize labels to -1
Next = zeros (nlma, 1) ;
Prev = zeros (nlma, 1) ;
nsets = 0 ;
nlabels = 0 ;
nedges = 0 ;

% pair up nodes that are close to each other
ave_dist = 0 ;
nave = 0 ;
for j = 1:nlma
    if ( Label (j) == -1 )
        R = LMA - LMA (:, j)*LMA_ones ;
        Rnorm = sum (R.^2) ;
        Rnorm (j) = inf ;
        [rmin k] = min (Rnorm) ;
        rmin = sqrt (rmin) ;
        if ( Label (k) == -1 ) % start a new set
            nsets = nsets + 1 ;
            Sets (nsets) = j ;
            Next (j) = k ;
            Prev (j) = k ;
            Next (k) = j ;
            Prev (k) = j ;
            Label (j) = nsets ;
            Label (k) = nsets ;
        else                    % add j to set Label (k) ;
            l = Prev (k) ;
            Prev (j) = l ;
            Next (l) = j ;
            Next (j) = k ;
            Prev (k) = j ;
            Label (j) = Label (k) ;
        end
        nedges = nedges + 1 ;
        EdgeList (1, nedges) = j ;
        EdgeList (2, nedges) = k ;
    end
end
if ( ~Quiet )
    fprintf (1, 'Initial number of edges: %i\n', nedges) ;
    fprintf (1, 'Initial number of sets: %i\n', nsets) ;
end

ActiveSet = ones (nsets, 1) ;
nactive = nsets ;
AliveNodes = ones (nlma, 1) ;
while ( nactive > 1 ) % add edges while no sets are isolated
%   fprintf (1, 'begin loop, nactive: %i\n', nactive) ;
    i = 0 ;
    while ( i <= nsets )
        length(ActiveSet);
        nactive;
        i = i + 1 ;
%       fprintf (1, 'set # %i nsets: %i\n', i, nsets) ;
        if ( ActiveSet (i) ) % the set is not isolated from the rest
%           fprintf (1, '    active\n') ;
            j = Sets (i) ;
            jstart = j ;
            LMA_ones = AliveNodes ;
            while ( 1 )
                j = Next (j) ;
                LMA_ones (j) = 0 ;
                if ( j == jstart )
                    break ;
                end
            end
            Comp = find (LMA_ones) ; % find the nodes in the complement
            lenComp = length (Comp) ;
%           fprintf (1, 'size of complement: %i\n', lenComp) ;
            LMAComp = LMA (:, Comp) ;
            BestDist = inf ;
            ones_Comp = ones (1, lenComp) ;
            while ( 1 )
                j = Next (j) ;
                R = LMAComp - LMA (:, j)*ones_Comp ;
                Rnorm = sqrt (sum (R.^2)) ;
                [rmin k] = min (Rnorm) ;
                if ( rmin < BestDist )
                    BestDist = rmin ;
                    BestK = Comp (k) ;
                    BestJ = j ;
                end
                if ( j == jstart )
                    break ;
                end
            end
%           fprintf (1, 'dist to adjacent set: %e\n', BestDist) ;
            if ( BestDist <= channel_connect ) % connect two channel segments
                nactive = nactive - 1 ;
                nedges = nedges + 1 ;
                EdgeList (1, nedges) = BestJ ;
                EdgeList (2, nedges) = BestK ;
                l = Label (BestK) ;
                if ( l < nsets )    % relabel the last set as l
                    m = Sets (nsets) ;
                    Sets (l) = m ;
                    ActiveSet (l) = ActiveSet (nsets) ;
                    n = m ;
                    while ( 1 )
                        n = Next (n) ;
                        Label (n) = l ;
                        if ( n == m )
                            break ;
                        end
                    end
                end
                nsets = nsets - 1 ;

                % update the labels
                j = BestK ;
                while ( 1 )
                    j = Next (j) ;
                    Label (j) = Label (BestJ) ;
                    if ( j == BestK )
                        break ;
                    end
                end

                % update the links
                m = Prev (BestK) ;
                n = Next (BestJ) ;
                Next (BestJ) = BestK ;
                Prev (BestK) = BestJ ;
                Next (m) = n ;
                Prev (n) = m ;
            else                    % this set is separated from the rest
%               fprintf (1, '    kill\n') ;
                ActiveSet (i) = 0 ;  % this set is dead
                nactive = nactive - 1 ;
                j = jstart ;
                while ( 1 )
                    j = Next (j) ;
                    AliveNodes (j) = 0 ;
                    if ( j == jstart )
                        break ;
                    end
                end
            end
        end
    end
end

% Remove dangling nodes and their associated edge. Changed by J.Lapierre

if ( RemoveDangle )
    degree = zeros (nlma, 1) ;
    adj = zeros (nlma, 1) ;
    edge = zeros (nlma, 1) ;
    for j = 1:nedges
        k = EdgeList (1, j) ;
        l = EdgeList (2, j) ;
        adj (k) = l ;
        adj (l) = k ;
        edge (k) = j ;
        edge (l) = j ;
        degree (k) = degree (k) + 1 ;
        degree (l) = degree (l) + 1 ;
    end
    J = find (degree == 1) ;
    EdgeFlag = ones (nedges, 1) ;
    ndangling = 0 ;
    for j = J'
        k = adj (j) ;
        if ( degree (k) > 2 )
            degree (k) = degree (k) - 1 ;
            EdgeFlag (edge (j)) = 0 ;
            ndangling = ndangling + 1 ;
        end
    end
    if ( ~Quiet )
        fprintf (1, 'number of dangling nodes: %i\n', ndangling) ;
    end
%     if ( RemoveDangle )
        K = find (EdgeFlag) ;
        EdgeList = EdgeList (:, K) ;
        nedges = size (EdgeList, 2) ;

        % remove dangling sets and their associated edge
        ndrop = 0 ;
        DropList = zeros (nsets, 1) ;
        for i = 1:nsets
            j = Sets (i) ;
            jstart = j ;
            n = 0 ;
            while ( 1 )
                j = Next (j) ;
                n = n + 1 ;
                if ( j == jstart )
                    break ;
                end
            end
            if ( n < DanglingChannel ) % drop the set
                ndrop = ndrop + 1 ;
                DropList (ndrop) = i ;
            end
        end
        if ( ~Quiet )
            fprintf (1, 'number of dangling sets: %i\n', ndrop) ;
        end
        EdgeLabels = zeros (2, nedges) ;
        for k = 1:nedges
            EdgeLabel (1, k) = Label (EdgeList (1, k)) ;
            EdgeLabel (2, k) = Label (EdgeList (2, k)) ;
        end
        EdgeFlag = ones (nedges, 1) ;
        for i = 1:ndrop
            j = DropList (i) ; % drop this set
            J = find (EdgeLabel (1, :) == j) ;
            EdgeFlag (J) = 0 ;
            J = find (EdgeLabel (2, :) == j) ;
            EdgeFlag (J) = 0 ;
        end
        K = find (EdgeFlag) ;
        EdgeList = EdgeList (:, K) ;
        nedges = size (EdgeList, 2) ;
%     end
else
    EdgeList = EdgeList (:, 1:nedges) ;
end

if ( ~Quiet )
    fprintf (1, 'final number of edges: %i\n', nedges) ;
end

if ( nargout > 1 )
    Vertex1 = zeros (3, nedges) ;
    Vertex2 = zeros (3, nedges) ;
    for j = 1:nedges
        Vertex1 (:, j) = LMA (:, EdgeList (1, j)) ;
        Vertex2 (:, j) = LMA (:, EdgeList (2, j)) ;
    end
end

if ( ~PlotGraph && ~Rotate )
    % map back to original nodes
    for j = 1:nedges
        EdgeList (1, j) = LMA_keep (EdgeList (1, j)) ;
        EdgeList (2, j) = LMA_keep (EdgeList (2, j)) ;
    end
    return ;
end

if ( PlotGraph && ~Rotate )
    Camera_frames = 1 ;
end

dlon = (lon_east - lon_west)/map_dimx;   %Degrees longitude per topo cell
dlat = (lat_north - lat_south)/map_dimy; %Degrees latitude per topo cell
load (loadmap) ;
topo = reshape (topo, map_dimx, map_dimy) ;
maxLat = max (LMAlat) ;
minLat = min (LMAlat) ;
maxLon = max (LMAlon) ;
minLon = min (LMAlon) ;
maxLat = maxLat + TopoExpandY/Flat ;    % max Lat in the plot
minLat = minLat - TopoExpandY/Flat ;    % min Lat in the plot
maxLon = maxLon + TopoExpandX/Flon ;    % max Lon in the plot
minLon = minLon - TopoExpandX/Flon ;    % min Lon in the plot


% % Find the cells in the topographic map that are beneath the flash.
cell_north = floor ((lat_north - maxLat)/dlat) ;
cell_south =  ceil ((lat_north - minLat)/dlat) ;
cell_west = floor ((minLon - lon_west)/dlon) ;
cell_east = ceil ((maxLon - lon_west)/dlon) ;

% % Get altitude above each cell, check that the flash is inside the map
FlashInside = 1 ;
if ( (cell_north < 1) || (cell_south > map_dimy) || ...
      (cell_west < 1) || (cell_east > map_dimx) )
    fprintf (1, 'The flash extends outside the specified topographic map\n')
    FlashInside = 0 ;
    west_cutoff = max([cell_west 1]) ;
    east_cutoff = min([cell_east map_dimx]) ;
    north_cutoff = max([cell_north 1]) ;
    south_cutoff = min([cell_south map_dimy]) ;
    Zgrid = topo (west_cutoff:east_cutoff, north_cutoff:south_cutoff)'/1000 ;%km
else
    Zgrid = topo (cell_west:cell_east, cell_north:cell_south)'/1000 ; % km
end

clear topo ;

% % Find the X and Y coordinates of the cells relative to the specified POI.
X = (lon_west - POI_lon + [cell_west:cell_east]*dlon)*Flon/1000 ;    % km
Y = (lat_north - POI_lat - [cell_north:cell_south]*dlat)*Flat/1000 ; % km
% if ( FlashInside )
%     nY = length (Y) ;
%     Xgrid = ones (nY, 1)*X ;
%     nX = length (X) ;
%     Ygrid = Y'*ones (1, nX) ;
% else
%     Xcutoff = (lon_west - POI_lon + ...
%               [west_cutoff:east_cutoff]*dlon)*Flon/1000 ;    % km
%     Ycutoff = (lat_north - POI_lat - ...
%               [north_cutoff:south_cutoff]*dlat)*Flat/1000 ; % km
%     nY = length (Ycutoff) ;
%     Xgrid = ones (nY, 1)*Xcutoff ;
%     nX = length (Xcutoff) ;
%     Ygrid = Ycutoff'*ones (1, nX) ;
% end
% 
figure(FigureHandle) ;
clf ;
% if ( NoAxes )
%     axis off
% end
% if ( (length (Xgrid) > 0) && (length (Ygrid) > 0) )
% %     surf(Xgrid,Ygrid,Zgrid)
% %     return
% %     mapshow (Xgrid, Ygrid, Zgrid)%, 'DisplayType', 'surface') ;
%     if ( NoAxes )
%         axis off
%     end
%     view(3)
% %     colormap(Zgrid)
%     grid on
% end
hold on

%LMA data relative to POI
LMAx = (LMAlon-POI_lon)*Flon/1000 ; % km
LMAy = (LMAlat-POI_lat)*Flat/1000 ; % km
LMAz = LMAalt/1000 ;                % km

Xmin = min(X) ;
Xmax = max(X) ;
Ymin = min(Y) ;
Ymax = max(Y) ;
Xmid = .5*(Xmin + Xmax) ;
Ymid = .5*(Ymin + Ymax) ;
Zmid = .5*(min(LMAalt(:)) + max(LMAalt(:)))/1000 ;

% Set up the camera direction and position
r = norm ([Xmax-Xmin Ymax-Ymin]) ; % length of diagonal
Camera_radius = r + Camera_dr/1000 ;       % camera radius
camtarget([Xmid, Ymid, Zmid]);             % camera target
camva(35) ;                                % camera viewing angle
cmap = load ('AnimateCmap') ;
ncolors = size (cmap, 1) ;
Camera_alt = (max(LMAz) + Camera_dz/1000) ; % altitude of view in km
% 
% Zmin = min (LMAz) ;
% if ( (length (Xgrid) > 0) && (length (Ygrid) > 0) )
%     tmin = min([Zmin min(min(Zgrid))]) ;
% end
% 
% xlim ([Xmin Xmax]) ;
% ylim ([Ymin Ymax]) ;
% zlim ([tmin max(LMAz)]) ;
% 
% % If no axes, then mark 2 km line segments on the southwest corner of the plot.
% if ( NoAxes )
%     x = Xmin ;
%     y = Ymin ;
%     if ( (length(Xgrid) > 0) && (length(Ygrid) > 0) )
%         z = Zgrid (end, 1) ;
%     else
%         z = Zmin ;
%     end
%     plot3 ([x x+2], [y y], [z+.1 z+.1], '-k', 'LineWidth', 3) ;
%     plot3 ([x x], [y y+2], [z+.1 z+.1], '-k', 'LineWidth', 3) ;
%     plot3 ([x x], [y y], [z+.1 z+2.1], '-k', 'LineWidth', 3) ;
%     text (x, y, z-.55, '2 km') ;
%     text (x, y+2.5, z+.1, 'N') ;
%     text (x+2.5, y, z+.1, 'E') ;
% else
%     if ( (length(Xgrid) > 0) && (length(Ygrid) > 0) )
%         x = Xmin ;
%         y = Ymax ;
%         z = Zgrid (1, 1) ;
%         text (x, y, z+1, 'N', 'FontSize', 16) ;
%         x = Xmax ;
%         y = Ymin ;
%         z = Zgrid (end, end) ;
%         text (x, y, z+1, 'E', 'FontSize', 16) ;
%     end
% end

% Plot the pulse graph
LMA = LMA/1000 ;
for j = 1:nedges
    k = EdgeList (1, j) ;
    l = EdgeList (2, j) ;
    x1 = LMA (1, k) ;
    x2 = LMA (1, l) ;
    X = [ x1 x2 ] ;
    y1 = LMA (2, k) ;
    y2 = LMA (2, l) ;
    Y = [ y1 y2 ] ;
    z1 = LMA (3, k) ;
    z2 = LMA (3, l) ;
    Z = [ z1 z2 ] ;
    plot3 (X, Y, Z, '-k', 'LineWidth', 2) ;
end

x = Camera_radius*cos(Camera_initial_angle) + Xmid ;
y = Camera_radius*sin(Camera_initial_angle) + Ymid ;
campos([x, y, Camera_alt]) ;
drawnow
% pause(3) ;

% % Plot the LMA points
% dt = (LMAt(end)-LMAt(1))/ncolors ;
% LMAc = 1 + floor ((LMAt - LMAt(1))/dt) ;
% I = find (LMAc > ncolors) ;
% LMAc (I) = ncolors ;
% LMAc = cmap(LMAc, :) ;
% nlma = length (LMAx) ;
% for i = 1:nlma
%     plot3 (LMAx(i), LMAy(i), LMAz(i), 'o', ...
%     'Markersize', 6, 'color', LMAc(i, :)) ;
% end
% drawnow
% % for unknown reasons, scatter3 changes the colors of the earth
% %scatter3 (LMAx, LMAy, LMAz, 11, LMAc, 'filled') ;

if ( ~Rotate )
    return
end
axis vis3d
theta = 2*pi/Camera_frames ;
frame_number = 0 ;
for n = 1:Camera_frames
    x = Camera_radius*cos(Camera_initial_angle + theta*n) + Xmid ;
    y = Camera_radius*sin(Camera_initial_angle + theta*n) + Ymid ;
    campos([x, y, Camera_alt]) ;
    drawnow
%   M(frame_number) = getframe ;
    if ( Avi )
        frame_number = frame_number + 1 ;
        M(frame_number) = getframe(gcf) ;
    end
end
if ( Avi )
    if ( length (AviFilename) > 4 )
        AviFilename
        if ( strcmp ('.avi', AviFilename (end-3:end)) )
            movie2avi (M, AviFilename, 'fps', Camera_fps) ;
        else
            movie2avi (M, strcat (AviFilename, '.avi'), 'fps', Camera_fps) ;
        end
    else
        movie2avi (M, strcat (AviFilename, '.avi'), 'fps', Camera_fps) ;
    end
end
