function LCS(solver, PEFile, startDate_start, startDate_end, startDay_shift, timeWindow, whichFTLE, plot_depth_idx, clim, compute_traj)

% Changable Inputs Begin
% ----------------------

accuracy = 'low';
interp_density = 1;

DataReadScript = 'read_aux_vars_from_data_pe_flat_NSF_Alpha';

SaveFig = 1;
SaveData = 1;
clip_axis = 1;

PlotScripts{1} = 'PlotScript_alpha_t';
input_filename = 'start.txt';

CFL = 0.8;
dt = 60; % edit and make larger 

dateFormat = 'dd_mmmm_HH_MM_SS';


% Changable Inputs End
% --------------------

if nargin == 0
    disp('Running with default values');
    solver = 1;
    PEFile = '/home/phaley/Projects/MassBay/2019/Jul15/Extract/pe_out_2_140m_MBAY_w_long.nc';
    startDate_start = [2001 6 7 0 0 0];
    startDate_end = [2001 6 7 0 0 0];
    startDay_shift = Inf;
    timeWindow = 19;
    whichFTLE = 'backward';
    plot_depth_idx = [];
    clim = [];
end
if nargin == 9
    compute_traj = 0;
end
warning off;

% Generate and add required paths
% -------------------------------

parentDir_id = 'LCS';
workDir_id = 'WorkDir';
dataDir_id = 'Data';
caseDir_id = 'Cases';

currDir = pwd;
parentPath = currDir(1:(strfind(currDir, parentDir_id) + length(parentDir_id) - 1));
parentDataDir = strcat(currDir(1:(strfind(currDir, parentDir_id) + length(parentDir_id) - 1)), '/', dataDir_id, '/');
GroupRunName = currDir(strfind(currDir, caseDir_id) + length(caseDir_id) + 1:end);

addpath(genpath(sprintf('%s/Src', parentPath)));

% Create (if required) and read from the mat file of the simulation
% -----------------------------------------------------------------

pefile_name_cell = strsplit(PEFile, '/');
DataDir = parentDataDir(1:end-1);

maskfile_name = strcat('aux_', pefile_name_cell{end-2}, '_', pefile_name_cell{end-1});
if solver == 1 && interp_density ~=1
    maskfile_name = strcat(maskfile_name, '_interp_', num2str(interp_density));
end

if ~exist(DataDir, 'dir')
    mkdir(DataDir);
end

disp('Creating Mask from PE File');

if ~exist(sprintf('%s/%s.mat', DataDir, maskfile_name), 'file')
    fprintf('Creating MATLAB Mask File for This PE File\n\n');
    eval(['[app] = ', DataReadScript, '(PEFile, 0, 0);']);
    if solver == 1 && interp_density ~= 1
        app = interpolate_fields(app, interp_density);
    end
    %save(sprintf('%s/%s', DataDir, maskfile_name), 'app');
else
    fprintf('MATLAB Mask File Already Exists for This PE File\n\n');
    load(sprintf('%s/%s.mat', DataDir, maskfile_name));
end

app_pe = app;

% Compute start day and end day offsets
% -------------------------------------

startDate_start(3) = startDate_start(3) + floor(startDate_start(4)/24);
startDate_start(4) = mod(startDate_start(4), 24);

startDate_end(3) = startDate_end(3) + floor(startDate_end(4)/24);
startDate_end(4) = mod(startDate_end(4), 24);

fprintf('Start Date of the PE Simulation is %s\n\n', datestr(app.start_date));

if datenum(startDate_start) < datenum(app.start_date) || datenum(startDate_end) + timeWindow > datenum(app.start_date) + app.time(end)/86400
    fprintf('Start and end dates of the PE file are %s and %s\n', datestr(app.start_date, 0), datestr(datenum(app.start_date) + app.time(end)/86400, 0));
    fprintf('Your requested dates for the simulation are %s and %s\n', datestr(startDate_start, 0), datestr(datenum(startDate_end) + timeWindow, 0));
    error('The Requested Dates Are Not Contained in the PE File');
else
    startDay_start = datenum(startDate_start) - datenum(app.start_date);
    startDay_end = datenum(startDate_end) - datenum(app.start_date);
end

% Setup
% -----

app.GroupRunName = GroupRunName;

app.ncid = netcdf(app.ncFile);

app.accuracy = accuracy;
app.interp_density = interp_density;

app.SaveFig = SaveFig;
app.SaveData = SaveData;
app.clip_axis = clip_axis;
app.clim = clim;

app.timeWindow = timeWindow;
app.plot_depth_idx = plot_depth_idx;

app.PlotScripts = PlotScripts;

app.CFL = CFL;

app.type = 1;
app.dim = 2;

app.whichFTLE = whichFTLE;

app.t0_sim = app.t0;
app.tEnd_sim = app.tEnd;
app.dt_sim = app_pe.dt;

if exist('dt', 'var')
    app.dt = dt;
else
    app.dt = 0.2*app.CFL/app.interp_density*min(app_pe.dx, app_pe.dy);	% Assume the velocity does not exceed 5 m/s
end


%% Save directory setup

app.ParentSaveDir = sprintf('%s/%s/%s', parentPath, workDir_id, app.GroupRunName);
app.ParentDir = app.ParentSaveDir;
app.SaveDirMat = sprintf('%s/mat', app.ParentSaveDir);
app.SaveDirImg = sprintf('%s/img', app.ParentSaveDir);

if ~exist(app.ParentSaveDir, 'dir')
    mkdir(app.SaveDirMat);
    mkdir(app.SaveDirImg);
end



%% Select solver, create grids
% ---------------------------

fprintf('Using a time step of %0.3f\n\n', app.dt);

if solver == 1
    fprintf('Using PDE based computation\n\n');
    app.solver = 'flowmap_2D';
elseif solver == 2
    fprintf('Using trajectory based computation\n\n');
    app.solver = 'flowmap_traj_3D';
end

eval('define_grids');

app_params = app;
lon = app.lon;
lat = app.lat;

%

% Start computation for all the required start / end dates
% --------------------------------------------------------

fprintf('Starting Flow Map and FTLE Computation\n\n');

for i = startDay_start : startDay_shift : startDay_end
    
    app = catstruct(app_pe, app_params);
    
    % Compute offsets, set up run name and (sub)directories
    % -----------------------------------------------------
    
    app.startDay = i;
    app.endDay = i + timeWindow;
    
    startDate_str = datestr(datenum(app.start_date) + app.startDay, dateFormat);
    endDate_str = datestr(datenum(app.start_date) + app.endDay, dateFormat);

    app.startDate_str = startDate_str;
    app.endDate_str = endDate_str;
    
    fprintf('Start Date is %s\n\n', startDate_str);
    
    app.RunName = sprintf('%s_Start_%s_End_%s', app.whichFTLE, startDate_str, endDate_str);
    app.SaveDir = sprintf('%s/%s', app.SaveDirMat, app.RunName);
    
    app.t0 = app.t0_sim + 86400*app.startDay;
    app.tEnd = app.t0_sim + 86400*app.endDay;
    
    if (~exist(app.SaveDir, 'dir'))
        mkdir(app.SaveDir);
    end
    
    % Compute and plot and save
    % -------------------------
    
    eval(['[app, final_flowmap] = ', app.solver, '(app);']);
    save(strcat(app.SaveDir, '/fmap.mat'), 'final_flowmap', 'lon', 'lat');
	init_conc_type=5;
    alpha_t = init_conc(app, init_conc_type, final_flowmap.x, final_flowmap.y); 
    %[pd, ftle, pra] = compute_ftle(app, final_flowmap);
    %save(strcat(app.SaveDir, '/ftle.mat'), 'final_flowmap', 'pd', 'ftle', 'lon', 'lat');

     for j = 1:length(app.PlotScripts)
       eval(app.PlotScripts{j});
     end
 
end

disp('All done!');
