function [geosap] = geosap(fname)
%GEOSAP Creates statistics and plots based on a CSV input from GEOSAP.
%
%   See also LLA2ENU.

fid = fopen(fname,'r');
inputs = textscan(fid,'%s%s%s%s%s%s%n%n%n%n%n%n','HeaderLines',1,'Delimiter',',');
status = fclose(fid); %#ok<NASGU>

geosap = create_structure(inputs);

states = fieldnames(geosap);
for istate = 1:length(states)
    areas = fieldnames(geosap.(states{istate}));
    for iarea = 1:length(areas)
        area = geosap.(states{istate}).(areas{iarea});
        % replace area with a version to which statistics have been added
        geosap.(states{istate}).(areas{iarea}) = process_area(area);
    end
end

[pathstr,name,ext,versn] = fileparts(fname); %#ok<NASGU,NASGU>
fid = fopen(fullfile(pathstr,[name '.txt']),'w');
states = fieldnames(geosap);
area_fields = fieldnames(geosap.(states{1}).(areas{1}));
remove_fields = {'lifts' 'trail_centerlines' 'trail_boundaries' 'boundary' 'lift_count' 'trail_count'};
statnames = setdiff(area_fields,remove_fields);
for istate = 1:length(states)
    areas = fieldnames(geosap.(states{istate}));
    for iarea = 1:length(areas)
        area = geosap.(states{istate}).(areas{iarea});
        
        count = fprintf(fid,'Area statistics for %s, %s:\n',areas{iarea},states{istate}); %#ok<NASGU>
        count = fprintf(fid,' lift_count = %g\n',area.lift_count);
        count = fprintf(fid,' trail_count = %g\n',area.trail_count);
        for istat = 1:length(statnames)
            statname = statnames{istat};
            value = m2ft(area.(statname));
            count = fprintf(fid,' %s = %g\n',statname,value); %#ok<NASGU>
        end
        
        count = fprintf(fid,'Individual lift lengths, vertical:\n'); %#ok<NASGU>
        for ilift = 1:area.lift_count
            n = area.lifts(ilift).name;
            l = m2ft(area.lifts(ilift).length);
            v = m2ft(area.lifts(ilift).vertical);
            count = fprintf(fid,' %s = %g, %g\n',n,l,v); %#ok<NASGU>
        end
        
        count = fprintf(fid,'Individual trail lengths, vertical:\n'); %#ok<NASGU>
        for itrail = 1:area.trail_count
            n = area.trail_centerlines(itrail).name;
            l = m2ft(area.trail_centerlines(itrail).length);
            v = m2ft(area.trail_centerlines(itrail).vertical);
            count = fprintf(fid,' %s = %g, %g\n',n,l,v); %#ok<NASGU>
        end
    end
end
status = fclose(fid); %#ok<NASGU>



function [geosap] = create_structure(inputs)
state = inputs{1};
area = inputs{2};
type = inputs{3};
group = inputs{4};
name = inputs{5};
desc = inputs{6};
color = inputs{7:9};
lla = [inputs{10:12}];
% convert lat/lon/alt to radians/radians/meters
%lla(:,1:2) = deg2rad(lla(:,1:2));       % requires mapping toolbox?
%lla(:,3) = distdim(lla(:,3),'ft','m');  % requires mapping toolbox?
lla(:,1:2) = pi./180.*lla(:,1:2);
lla(:,3) = ft2m(lla(:,3));

%% desired data format
% geosap.state().area().lifts().name
% geosap.state().area().lifts().description
% geosap.state().area().lifts().color
% geosap.state().area().lifts().lla
% geosap.state().area().trail_centerlines().name
% geosap.state().area().trail_centerlines().lla

% unique 'rows' option doesn't work on 2D cell arrays, so must concatenate
% state & area
[ustatearea,ind] = unique(cell2mat([state area]),'rows');
ustatearea = [state(ind) area(ind)];
% remove spaces so state & area can be used as fielnames
ustatearea = strrep(ustatearea,' ','_');
ind = [1 ind];

isl = strcmp(type,'Lifts');
istc = strcmp(type,'Trail Centerlines');
istb = strcmp(type,'Trail Boundaries');
isb = strcmp(type,'Boundary');

for iarea = 1:size(ustatearea,1)
    area_ind = ind(iarea):ind(iarea+1);
    area = struct('lifts',[],'trail_centerlines',[],'trail_boundaries',[],'boundary',[]);
    
    area.lifts = create_lifts(name(area_ind(isl)),lla(area_ind(isl),:));
    area.trail_centerlines = create_trail_centerlines(name(area_ind(istc)),lla(area_ind(istc),:));
%     area.trail_boundaries = create_trail_boundaries(name(area_ind(istb)),lla(area_ind(istb),:));
%     area.boundary = create_boundary(name(area_ind(isb)),lla(area_ind(isb),:));
    
    geosap.(ustatearea{iarea,1}).(ustatearea{iarea,2}) = area;
end



function [lifts] = create_lifts(name,lla)
% TODO: figure out if there's a way to avoid needing unique lift names
% while still allowing midstation/bend points - maybe a separating row?

% For now, assume that 2 lifts in a row will have unique names
% This allows duplicate named lifts where there is
% another lift in the middle.

% could do some calculations and see if points are too far apart to be one lift,
% but this is going to happen most often with short lifts (e.g. unnamed surface tows),
% for which this check would be ineffective.  could also check for lift
% crossings, but this is also not foolproof nor worth the complication.

% TODO: handle no points, 1 point cases

% write the first two points
lifts.name = name{1};
lifts.lla = lla(1:2,:);
ilift = 1;
ipoint = 3;
while ipoint <= length(name)
    if strcmp(name(ipoint),name(ipoint-1))
        lifts(ilift).lla = [lifts(ilift).lla; lla(ipoint,:)];
    else
        ilift = ilift + 1;
        lifts(ilift).name = name{ipoint};
        lifts(ilift).lla = lla(ipoint:(ipoint+1),:);
        ipoint = ipoint + 1;
    end
    ipoint = ipoint + 1;
end



function [trail_centerlines] = create_trail_centerlines(name,lla)
% TODO: unify processing that duplicates create_lifts
% write the first two points
trail_centerlines.name = name{1};
trail_centerlines.lla = lla(1:2,:);
itrail = 1;
ipoint = 3;
while ipoint <= length(name)
    if strcmp(name(ipoint),name(ipoint-1))
        trail_centerlines(itrail).lla = [trail_centerlines(itrail).lla; lla(ipoint,:)];
    else
        itrail = itrail + 1;
        trail_centerlines(itrail).name = name{ipoint};
        trail_centerlines(itrail).lla = lla(ipoint:(ipoint+1),:);
        ipoint = ipoint + 1;
    end
    ipoint = ipoint + 1;
end



% process area
function [area] = process_area(area)
% eventually this should be vectorized, but for now, add ENU to each lift

% find base point, defined as the lowest altitude lift point
[base_alt,ind] = min(area.lifts(1).lla(:,3));
base = area.lifts(1).lla(ind,:);
for ilift = 2:length(area.lifts)
    [lift_alt,ind] = min(area.lifts(ilift).lla(:,3));
    if lift_alt < base_alt
        base_alt = lift_alt;
        base = area.lifts(ilift).lla(ind,:);
    end
end

% for each lift, add east-north-up cartesian system about the base
for ilift = 1:length(area.lifts)
    lla = area.lifts(ilift).lla;
    % lla2enu requires longitude to be 0 < x < 2pi
    lla(:,2) = mod(lla(:,2),2.*pi);
    area.lifts(ilift).enu = lla2enu(lla,base);
end

% for each trail, add east-north-up cartesian system about the base
for itrail = 1:length(area.trail_centerlines)
    lla = area.trail_centerlines(itrail).lla;
    % lla2enu requires longitude to be 0 < x < 2pi
    lla(:,2) = mod(lla(:,2),2.*pi);
    area.trail_centerlines(itrail).enu = lla2enu(lla,base);
end


%% Individual Area Statistics/Plots


%% Lift Statistics
% lift count
lift_count = length(area.lifts);
area.lift_count = lift_count;
% total length
% total vertical drop
% area vertical drop
lift_length = zeros(lift_count,1);
lift_minalt = zeros(lift_count,1);
lift_maxalt = zeros(lift_count,1);
for ilift = 1:lift_count
    lift = area.lifts(ilift);
    
    % loop through all points to make sure lifts that have midstations
    % and/or bends are treated properly
    cum_length = 0;
    for ipoint = 1:(size(lift.enu,1)-1)
        dist = sqrt(sum((lift.enu(ipoint,:) - lift.enu(ipoint+1,:)).^2));
        cum_length = cum_length + dist;
    end
    area.lifts(ilift).length = cum_length;
    area.lifts(ilift).vertical = max(lift.lla(:,3)) - min(lift.lla(:,3));
    
    lift_length(ilift) = cum_length;
    lift_minalt(ilift) = min(lift.lla(:,3));
    lift_maxalt(ilift) = max(lift.lla(:,3));
end

area.lift_length = sum(lift_length);
lift_vertical = lift_maxalt - lift_minalt;
area.lift_vertical = sum(lift_vertical);
area.vertical = max(lift_maxalt) - min(lift_minalt);


% properties-based:
% average speed (fpm, m/s)
    % take total length of all lifts, divided by total time to ride each lift once
% beginner-deweighted speed (fpm, m/s)
    % simple version: exclude all lifts which serve <25% of the max lift vertical 
    % fancy version: take total length of all lifts which serve >75% beginner terrain,
    % divided by total time to ride those lifts once
% uphill capacity (people/hr)
% uphill vertical (people-ft/hr)

% advanced:
% expected lines

%% Lift Plots
% vertical vs. horizontal run
% elevation vs. horizontal run

%% Trail Centerline Statistics
trail_count = length(area.trail_centerlines);
area.trail_count = trail_count;
% length
% vertical
% TODO: unify processing that is duplicated for lifts
for itrail = 1:trail_count
    trail = area.trail_centerlines(itrail);
    
    % loop through all points
    cum_length = 0;
    for ipoint = 1:(size(trail.enu,1)-1)
        dist = sqrt(sum((trail.enu(ipoint,:) - trail.enu(ipoint+1,:)).^2));
        cum_length = cum_length + dist;
    end
    area.trail_centerlines(itrail).length = cum_length;
    area.trail_centerlines(itrail).vertical = max(trail.lla(:,3)) - min(trail.lla(:,3));
end

area.trail_length = sum([area.trail_centerlines.length]);
area.trail_vertical = sum([area.trail_centerlines.vertical]);

% mean, percentile(10,25,50-median,75,90), max pitch
    % percentile pitch helps characterize how sustained the pitch is
% undulation (vertical pitch change)
% curvature (horizontal direction change)
    % need to disregard the effect of minor zig zags (use percentage of
    % slope width to detect)

%% Trail Centerline Plots


%% Trail Boundary Statistics

%% Trail Boundary Plots

%% Area Boundary Statistics

%% Area Boundary Plots

%% Multi-Area Plots
% total vertical


%% Lift Statistics

%% Lift Plots

%% Trail Centerline Statistics

%% Trail Centerline Plots

%% Trail Boundary Statistics

%% Trail Boundary Plots

%% Area Boundary Statistics

%% Area Boundary Plots



%% Error Checking
function [] = validate_area()
% color validation
% trail centerline/boundary names match
% trail centerline/boundary order matches
% each lift has > 1 point
% each trail centerline has > 1 point
% each trail centerline does not cross itself
% each trail boundary has > 3 points
% each trail boundary does not cross itself
% boundary includes all others horizontally


%% Utilities
% Splice trails
% Trails served by lift


function [ft] = m2ft(m)
ft = m./0.3048;

function [m] = ft2m(ft)
m = ft.*0.3048;