function returnResult = PlotPhaseDiagrams(regimeStruct, times, vars, varBounds, DEqns_comp, ...
    allDEixMap, numExt, networkObjectDir, DSSRTfigHandle, verboseTog)
% Version 1.0, (c) Robert Clewley, Cornell University, 2005
%
% Plot sequence of phase-plane diagrams from reduced dynamical regimes
% (generated by RegimeDet), or produce movies

% -------------------------------------------------------------------------

%%% TEMP DEBUGGING %%%
verboseTog = true;
disp('TEMP DEBUGGING: `verboseTog` forced true')
disp(' ')

% -------------------------------------------------------------------------

% INFO on regime structure's fields, on entry to this function:
%regimeStruct.regimes
%regimeStruct.numEps
%regimeStruct.numRegs
%regimeStruct.varnames
%regimeStruct.transSeq
%regimeStruct.isCycle
%regimeStruct.period
%regimeStruct.startTime    % not necessarily the same as the value passed to RegimeDet
%regimeStruct.tScaleThresh
%regimeStruct.algOptions
%regimeStruct.focVarIx

% INFO on regime structure's fields:
% epochs, timeInt, dynVars, dynCross, nonDynVars, fastVars, slowVars, dimensionA, dimensionB,
%   dimensionC, qsPars

% -------------------------------------------------------------------------

returnResult = [];

sysStrRoot='PlotPhaseDiagrams:  ';
sysStrRootIndent = length(sysStrRoot);

if nargin == 9
    verboseTog = false;
elseif nargin ~= 10
    disp([sysStrRoot,'Wrong number of arguments passed. Expected 9 or 10.'])
    return
end

%% Load global constants
GlobalConstants % this line will be replaced by the contents of this script file when compiled (see DSSRT.m)

disp([sysStrRoot,'Calculation may take some time while nullclines are calculated'])
disp('              (press <escape> to cancel)')
if verboseTog
    disp(' ')
end

% -------------------------------------------------------------------------

%% USER-SPECIFIED CONSTANTS
result = DialogBox(['number of quasi-static frames per regime, for slow vars (> 1 and < 1000 in magnitude, negative vals selects movie)'],5,verboseTog);
if ~isempty(result)
    splitRes = floor(str2num(result));
    if splitRes < 0
        makemovie = true;
        splitRes = abs(splitRes);
    else
        makemovie = false;
    end
    if abs(splitRes) <= 1 && abs(splitRes) >= 1000
        splitRes = 5;
        disp(' Value out of range -- using default of 5')
    end
else
    if verboseTog
        disp('           Cancelled')
    end
    return
end

if makemovie
    zoomTraj = true;
    numNullPoints = 50;
    % need more magnification because short trajectory segment per frame
    extendFac_def = 5.0; % default
    result = DialogBox(['magnification factor (> 1.5 and < 1000)'],extendFac_def,verboseTog);
	if ~isempty(result)
        extendFac = str2num(result);
        if extendFac <= 1.5 && extendFac > 1000
            extendFac = extendFac_def;
            fprintf(' Value out of range -- using default of %.4f\n',extendFac)
        end
	else
        if verboseTog
            disp('           Cancelled')
        end
        return
	end
else
	result = DialogBox(['number of points on nullclines (> 10 and < 1000)'],50,verboseTog);
	if ~isempty(result)
        numNullPoints = str2num(result);
        if numNullPoints < 10 && numNullPoints > 1000
            numNullPoints = 50;
            disp(' Value out of range -- using default of 50')
        end
	else
        if verboseTog
            disp('           Cancelled')
        end
        return
	end
	if verboseTog
        disp([sysStrRoot 'Zoom to trajectory neighbourhood? (see dialog box)'])
	end
	ButtonName = questdlg('Zoom choice','Zoom to trajectory neighbourhood?','Yes','No','No');
	switch ButtonName
        case 'Yes'
            zoomTraj = true;
        case 'No'
            zoomTraj = false;
        otherwise % dialog cancel button
            disp('          Cancelled')
            return
	end
    extendFac = 2.0;
end
% for traj to occupy approx 1/extendFac of the phase plane in each
% direction, we halve the extendFac value because it's applied twice in
% each direction.
extendFac_half = extendFac / 2;

if verboseTog
    disp([sysStrRoot 'Automatically determine phase-plane variables? (see dialog box)'])
end
ButtonName = questdlg('Auto-choose variables','Automatically determine phase-plane variables','Auto','Manual','Auto');
switch ButtonName
    case 'Auto'
        autoAxes = true;
    case 'Manual'
        autoAxes = false;
    otherwise % dialog cancel button
        disp('          Cancelled')
        return
end

focusSet = regimeStruct.transSeq{TS_FOCUS};

if ~autoAxes
    defXix = regimeStruct.focVarIx;
    temp = setdiff(focusSet,defXix);
    defYix = temp(1);
	result = DialogBox(['x axis variable (defaults to focused var)'],regimeStruct.varnames{defXix},verboseTog);
	if ~isempty(result)
        Xname = result;
        [p pix] = ismember(Xname, regimeStruct.varnames);
        if pix==0
            Xname = regimeStruct.varnames{defXix};
            fprintf(' Invalid response -- using default of %s\n',Xname)
        end
        [p userXix] = ismember(Xname, regimeStruct.varnames);
	else
        if verboseTog
            disp('           Cancelled')
        end
        return
	end
    result = DialogBox(['y axis variable'],regimeStruct.varnames{defYix},verboseTog);
	if ~isempty(result)
        Yname = result;
        [p pix] = ismember(Yname, regimeStruct.varnames);
        if pix==0
            Yname = regimeStruct.varnames{defYix};
            fprintf(' Invalid response -- using default of %s\n',Yname)
        end
        [p userYix] = ismember(Yname, regimeStruct.varnames);
	else
        if verboseTog
            disp('           Cancelled')
        end
        return
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Phase 1: determine sequence of what will be drawn
% dimension based on dimensionB value
% if regime has > 2D look make fast and slow variables non-dynamic
% if regime still has > 2D ask which variables to plot in plane
% slow variables as quasi-static pars. -> sequence of plots for one regime, showing
%   phase-plane for equally-spaced values of the "parameter(s)"

% plotInfo is a flattened structure -- individual regimes and regimes with
% slowly varying parameters are flattened into one list of structures,
% having the following fields:
%Xix -- usually, the focused variable
%Xname
%Yix
%Yname
%regime
%timeInterval -- pair
%tposInterval -- pair
%traj -- list of pairs
%nullX -- list of pairs (full and reduced versions at t = start, halfway, end)
%nullY -- list of pairs (full and reduced versions at t = start, halfway, end)
%qsp -- cell array of cell arrays {varname, val}
%axesLimits -- list
%issplit -- boolean flag
%frame -- frame # for split regimes (0 otherwise)
%numFrames
%infoStr

plotInfo(1).qsp = []; % initial value
diagix = 1;
lastXix = 0;
lastYix = 0;
lastNumFrames = 0;
for regix = 1:regimeStruct.numRegs
    if makemovie
        disp([sysStrRoot 'Regime ' num2str(regix) ' - calculating sequence of nullclines and trajectory segments...'])
        drawnow
    end
    reg = regimeStruct.regimes(regix);
    doreg = true; % default
    dim = reg.dimensionB;
    if autoAxes
        Xix = regimeStruct.focVarIx;
        if dim == 1
            if reg.dimensionA > 1
                dim = reg.dimensionA;
                if dim > 2
                    % user selects Y axis variable (a slow or fast variable) ~= X
                    disp([sysStrRoot 'Regime ' num2str(regix) ': variables that can be viewed on Y-axis:'])
                    otherixs = union(reg.fastVars,reg.slowVars);
                    othernames = regimeStruct.varnames(otherixs);
                    for vix = reg.fastVars
                        fprintf('             %s (fast)\n',regimeStruct.varnames{vix});
                    end
                    for vix = reg.slowVars
                        fprintf('             %s (slow)\n',regimeStruct.varnames{vix});
                    end
                    result = DialogBox(['Y-axis variable (by name) for regime ' num2str(regix)], ...
                        othernames{1},verboseTog);
					if ~isempty(result)
                        Yname = result;
                        [p pix] = ismember(Yname, othernames);
                        if pix==0
                            Yname = othernames{1};
                            fprintf(' Invalid response -- using default of %s\n',Yname)
                        end
                        [p Yix] = ismember(Yname, regimeStruct.varnames);
                    else
                        return
					end
                else % dim == 2 (dim == 1 already covered in outer if statement)
                    Yix = setdiff(reg.dynVars,Xix);
                end
            else
                fprintf('     Regime %g has dimension 1 ("A"-measure) -- no diagram\n', regix)
                doreg = false;
            end
        elseif dim > 2
            % user selects Y axis variable ~= X
            getrid = [reg.fastVars(:)',reg.slowVars(:)'];
            othernames = regimeStruct.varnames(setdiff(union(reg.dynVars,reg.dynCross),getrid));
            disp([sysStrRoot 'Regime ' num2str(regix) ': variables that can be viewed on Y-axis:'])
            for vix = 1:dim-1
                fprintf('             %s\n', othernames{vix});
            end
            result = DialogBox(['Y-axis variable (by name) for regime ' num2str(regix)], ...
                othernames{1},verboseTog);
			if ~isempty(result)
                Yname = result;
                [p pix] = ismember(Yname, othernames);
                if pix==0
                    Yname = othernames{1};
                    fprintf(' Invalid response -- using default of %s\n',Yname)
                end
                [p Yix] = ismember(Yname, regimeStruct.varnames);
            else
                return
			end
        else % dim == 2, perfect!
            getrid = [Xix,reg.fastVars(:)',reg.slowVars(:)'];
            Yix = setdiff(union(reg.dynVars,reg.dynCross),getrid);
        end
    else
        Xix = userXix;
        Yix = userYix;
    end
    if ~doreg
        continue % regix for loop
    end
    % quasi-static "parameters" determine whether regime is split into
    % several diagrams
    if length(reg.qsPars) > 0
        split = splitRes;
    elseif length(setdiff(reg.slowVars,Yix)) > 0
        reg.qsPars = setdiff(reg.slowVars,Yix);  % for these purposes only
        split = splitRes;
    else
        split = 1;
    end

    % extract regime start and end time indices from the transition sequence assoc'd
    % with the regime
    numEpochs = length(reg.epochs);
    tposLo_reg = regimeStruct.transSeq{TS_TSEQ}{reg.epochs(1)}{TSEQ_POSN};
    if reg.epochs(numEpochs)+1 > length(regimeStruct.transSeq{TS_TSEQ})
        % get here generally for non-cycles
        tposHi_reg = regimeStruct.transSeq{TS_TSEQ}{1}{TSEQ_POSN} + regimeStruct.transSeq{TS_ENDP}-1;
    else
        % get here if TS for a cycle (when extra entry for end of period is
        % included)
        tposHi_reg = regimeStruct.transSeq{TS_TSEQ}{reg.epochs(numEpochs)+1}{TSEQ_POSN}-1;
    end
    currAbsActs = setdiff([reg.dynVars(:)',reg.dynCross(:)',reg.nonDynVars(:)',regimeStruct.focVarIx], ...
        [reg.fastVars(:)']);
    DEqns_comp = SetActivesSwitch(DEqns_comp, [Xix,Yix], allDEixMap, numExt, currAbsActs);
    
    % split regime time interval according to qsp's
    tstep = floor((tposHi_reg-tposLo_reg)/split);
    if tstep <= 1
        % usually for makemovie
        split_actual = tposHi_reg-tposLo_reg;
        tstep = 1;
    else
        split_actual = split;
    end
    
    tpos_halfway_reg = floor((tposHi_reg + tposLo_reg)/2);
    vars_halfway_reg = vars(tpos_halfway_reg,:); % record variable values at halfway point of regime
    
    for i = 1:split_actual
        % allow to refresh DSSRT and Matlab console windows periodically
        if mod(i,10) == 1
            drawnow
            if get(DSSRTfigHandle,'CurrentCharacter') == 27 % escape key
                beep
                disp([sysStrRoot 'Computation cancelled by user'])
                return
            end
        end
        tposLo = tposLo_reg+(i-1)*tstep;
        tposHi = min(tposLo_reg+i*tstep,tposHi_reg);
        tposLen = tposHi-tposLo;
        if tposLen < 1
            % then split resolution doesn't find us enough time points
            % so skip ahead
            if verboseTog
                fprintf('     Not enought time points. Skipping ahead at i = %g/%g\n', i, split_actual)
            end
            continue
        end
        
        % prepare plotInfo
        plotInfo(diagix).Xix = Xix;
        plotInfo(diagix).Xname = regimeStruct.varnames{Xix};
        plotInfo(diagix).Yix = Yix;
        plotInfo(diagix).Yname = regimeStruct.varnames{Yix};
        plotInfo(diagix).regime = regix;
        if split_actual > 1 || (lastXix==Xix && lastYix==Yix)
            if lastXix==Xix && lastYix==Yix
                lastFrame = lastNumFrames;
            else
                lastFrame = 0;
            end
            plotInfo(diagix).issplit = true;
            plotInfo(diagix).frame = i + lastFrame;
        else
            plotInfo(diagix).issplit = false;
            plotInfo(diagix).frame = i;
            lastFrame = 0;
        end
        lastXix = Xix;
        lastYix = Yix;
        
        % extract trajectory and nullclines
        DE_X = allDEixMap(Xix);
        DE_Y = allDEixMap(Yix);
        traj = [vars(tposLo:tposHi,Xix), vars(tposLo:tposHi,Yix)];
        if makemovie
            if exist('traj_temp','var') == 0
                traj_temp = [vars(tposLo_reg:tposHi_reg,Xix), vars(tposLo_reg:tposHi_reg,Yix)];
                maxtraj = max(traj_temp);
                mintraj = min(traj_temp);
            end % otherwise keep same values
        else
            maxtraj = max(traj);
            mintraj = min(traj);
        end
        extents = [abs(maxtraj(1)-mintraj(1)), abs(maxtraj(2)-mintraj(2))];
        boundsX = varBounds(Xix,:);
        boundsY = varBounds(Yix,:);
        % if zoomTraj option selected, add 100% to extent intervals for axis limits,
        % unless those exceed variable bounds
        if zoomTraj
            upperX = min(extents(1)*extendFac_half + mintraj(1), boundsX(2));
            lowerX = max(maxtraj(1) - extents(1)*extendFac_half, boundsX(1));
            upperY = min(extents(2)*extendFac_half + mintraj(2), boundsY(2));
            lowerY = max(maxtraj(2) - extents(2)*extendFac_half, boundsY(1));
        else
            upperX = boundsX(2);
            lowerX = boundsX(1);
            upperY = boundsY(2);
            lowerY = boundsY(1);
        end
        extents_new = [abs(upperX-lowerX), abs(upperY-lowerY)];
        
        % only prepare reduced system's nullclines if not for movie (in which case these
        % matrices will remain empty)
        % prepare Y nullcline
        tpos_halfway = floor((tposHi + tposLo)/2);
        varVals = vars(tpos_halfway,:);
        nullY_f = [];
        nullY_r = [];
        stepX = extents_new(1)/numNullPoints;
        for xval=lowerX:stepX:upperX
			varVals_est = PrepareFocVars(reg.fastVars, setdiff(reg.slowVars,reg.qsPars), varVals, vars_halfway_reg, ...
                Xix, xval, focusSet, DEqns_comp, allDEixMap);
            nyf = GetQsfpVal(DEqns_comp{allDEixMap(Yix)},varVals_est.full,false);
            if nyf < LARGEBOUND
                nullY_f = [nullY_f; [xval, nyf]];
            end
%             if ~makemovie
                nyr = GetQsfpVal(DEqns_comp{allDEixMap(Yix)},varVals_est.reduced,true);
                if nyr < LARGEBOUND
                    nullY_r = [nullY_r; [xval, nyr]];
                end
%             end
        end
        % prepare X nullcline
        varVals = vars(tpos_halfway,:);
        nullX_f = [];
        nullX_r = [];
        stepY = extents_new(2)/numNullPoints;
        for yval=lowerY:stepY:upperY
			varVals_est = PrepareFocVars(reg.fastVars, setdiff(reg.slowVars,reg.qsPars), varVals, vars_halfway_reg, ...
                Yix, yval, focusSet, DEqns_comp, allDEixMap);
            nxf = GetQsfpVal(DEqns_comp{allDEixMap(Xix)},varVals_est.full,false);
            if nxf < LARGEBOUND
                nullX_f = [nullX_f; [nxf, yval]];
            end
%             if ~makemovie
                nxr = GetQsfpVal(DEqns_comp{allDEixMap(Xix)},varVals_est.reduced,true);
                if nxr < LARGEBOUND
                    nullX_r = [nullX_r; [nxr, yval]];
                end
%             end
        end
        % fill in remaining fields of plotInfo
        plotInfo(diagix).timeInterval = [times(tposLo),times(tposHi)];
        if makemovie
            % start trajectory at beginning of regime for movies
            plotInfo(diagix).tposInterval = [tposLo_reg,tposHi];
            plotInfo(diagix).traj = traj_temp;
        else
            plotInfo(diagix).tposInterval = [tposLo,tposHi];
            plotInfo(diagix).traj = traj;
        end
        plotInfo(diagix).nullX_f = nullX_f;
        plotInfo(diagix).nullY_f = nullY_f;
        plotInfo(diagix).nullX_r = nullX_r;
        plotInfo(diagix).nullY_r = nullY_r;
        plotInfo(diagix).axesLimits = [lowerX upperX lowerY upperY];
        plotInfo(diagix).numFrames = split_actual + lastFrame;
        plotInfo(diagix).infoStr = '';
        plotInfo(diagix).qsp = {};
        if length(reg.qsPars) > 0
            for i = 1:length(reg.qsPars)
                if ~ismember(reg.qsPars(i), [Xix,Yix])
                    % e.g. if Xix and Yix were chosen by user, don't want
                    % these qsPars to clash with those choices
                    plotInfo(diagix).qsp = {plotInfo(diagix).qsp{:}, ...
                            {regimeStruct.varnames{reg.qsPars(i)}, varVals(reg.qsPars(i))}};
                end
            end
        end
        % continue i for loop
        lastNumFrames = plotInfo(diagix).numFrames; % for use in next regime, to combine regime movies
        diagix = diagix + 1;
    end
end

numdiags = length(plotInfo);

% -------------------------------------------------------------------------

%% Phase 2: Plot sequence of diagrams and info to console
% create figure w/ "phase diagram" + number/total for window title, and fig title
%   including qsp val and time interval
% plot IC and endpoint as crosses on trajectory
% plot trajectory segment
% Info to console:
%   number of diagrams appearing and list of info about
%     each (inc. title, time interval, qsp vals)

seqOffset = 20;
figTitle1 = 'Phase Diagram ';
date_string = strrep(datestr(now,1),'-','');

disp(' ')
fprintf( '   Number of diagrams: %g\n', numdiags)
if verboseTog
	disp('     Remember that the validity of the nullclines diminishes with distance')
	disp('     from the plotted trajectory segments, as different reduced dynamical')
	disp('     regimes are likely to apply far from the known orbit. The nullclines')
    disp('     are only approximate even very near the trajectory segment, as they')
    disp('     are calculated at an x-variable value half-way through the time interval')
    disp('     indicated on the associated diagram.')
	disp(' ')
end
disp(    '     (Full system`s nullclines = dashed, reduced system`s nullclines = dotted')
disp(    '     (x-nullcline = magenta, y-nullcline = blue)')

% graphics stuff
screenSize = get(0,'ScreenSize');
yheight = 375;
xwidth = 450;

if makemovie
    diagseq = 1:numdiags;
else
    diagseq = numdiags:-1:1;
    if verboseTog
        disp([sysStrRoot 'Displaying diagrams and their info in reverse order'])
    end
end

for diagix = diagseq
    % allow to refresh DSSRT and Matlab console windows periodically
    if mod(diagix,10) == 1
        drawnow
        if get(DSSRTfigHandle,'CurrentCharacter') == 27 % escape key
            beep
            disp([sysStrRoot 'Computation cancelled by user'])
            return
        end
    end
    % figure basics
    figTitle2 = [num2str(diagix) '/' num2str(numdiags) ' -> Regime ' num2str(plotInfo(diagix).regime)];
    if makemovie
        xbotleft = 50;
        ybotleft = screenSize(4)-yheight*1.2;
    else
        xbotleft = 50+seqOffset*(diagix-1);
        ybotleft = ceil(screenSize(4)-yheight*1.2-seqOffset*(diagix-1));
        if ybotleft < seqOffset % near bottom of screen
            xbotleft = xbotleft + seqOffset;
            doffset = diagix;
            ybotleft = screenSize(4)-yheight*1.2-seqOffset*(diagix-doffset);
        end
        if xbotleft > screenSize(3)-xwidth/2 % near right edge of screen
            xbotleft = ceil(screenSize(3)-xwidth/2);
        end
    end
    fighandle(diagix) = figure('NumberTitle', 'off', 'Name', [ figTitle1 figTitle2 ],...
            'Position', [xbotleft, ybotleft, xwidth, yheight]);
    tLo = plotInfo(diagix).timeInterval(1);
    tHi = plotInfo(diagix).timeInterval(2);
    tposLo = plotInfo(diagix).tposInterval(1);
    tposHi = plotInfo(diagix).tposInterval(2);
    timeIntervalStr = ['[',num2str(tLo),',',num2str(tHi),']'];
    if ~isempty(plotInfo(diagix).qsp)
        qspStr = 'Quasi-static:';
        for qix = 1:length(plotInfo(diagix).qsp)
            q = plotInfo(diagix).qsp{qix};
            if qix == length(plotInfo(diagix).qsp)
                commaStr = '';
            else
                commaStr = ',';
            end
            qspStr = [qspStr ' ' q{1} '=' num2str(q{2}) commaStr];
        end
    else
        qspStr = '';
    end
    infostr = ['t \in ' timeIntervalStr ' ' qspStr];
    title(infostr)
    axes(gca)
    set(gca,'position', [.15 .1 .8 .8], 'xlim', plotInfo(diagix).axesLimits(1:2), ...
        'ylim', plotInfo(diagix).axesLimits(3:4), 'linewidth', 1.0)
    xlabel(plotInfo(diagix).Xname)
    ylabel(plotInfo(diagix).Yname)
    hold on
    
    % print info to console
    if verboseTog
        if ~plotInfo(diagix).issplit || plotInfo(diagix).frame == 1
            % diagram for new regime
            if makemovie
                fprintf('\n   Making movie for regime %g\n', plotInfo(diagix).regime)
            else
                fprintf('\n   Diagram %g of %g -> regime %g:\n', diagix, numdiags, plotInfo(diagix).regime);
            end
            fprintf('     Axes: %s vs. %s\n', plotInfo(diagix).Xname, plotInfo(diagix).Yname);
            if ~makemovie
                disp([  '     ', infostr])
            end
            plotInfo(diagix).infoStr = infostr;
        end
    end

    % plot nullclines
    if ~isempty(plotInfo(diagix).nullX_f)
        plot(plotInfo(diagix).nullX_f(:,1),plotInfo(diagix).nullX_f(:,2),'m--','linewidth',2)
    elseif ~plotInfo(diagix).issplit && verboseTog
        disp(['     Note: full system`s ' plotInfo(diagix).Xname ' nullcline does not exist'])
    end
    if ~isempty(plotInfo(diagix).nullY_f)
        plot(plotInfo(diagix).nullY_f(:,1),plotInfo(diagix).nullY_f(:,2),'b--','linewidth',2)
    elseif ~plotInfo(diagix).issplit && verboseTog
        disp(['     Note: full system`s ' plotInfo(diagix).Yname ' nullcline does not exist'])
    end
    % only shows reduced system's nullclines if not for movie (in which case these
    % matrices are empty)
    if ~isempty(plotInfo(diagix).nullX_r)
        plot(plotInfo(diagix).nullX_r(:,1),plotInfo(diagix).nullX_r(:,2),'m:','linewidth',1.5)
    elseif ~plotInfo(diagix).issplit && verboseTog
        disp(['     Note: reduced system`s ' plotInfo(diagix).Xname ' nullcline does not exist'])
	end
    if ~isempty(plotInfo(diagix).nullY_r)
        plot(plotInfo(diagix).nullY_r(:,1),plotInfo(diagix).nullY_r(:,2),'b:','linewidth',1.5)
    elseif ~plotInfo(diagix).issplit && verboseTog
        disp(['     Note: reduced system`s ' plotInfo(diagix).Yname ' nullcline does not exist'])
	end
    
    % plot trajectory segment, with IC and endpoint crosses if not for movie
    if makemovie
        plot(plotInfo(diagix).traj(1:tposHi-tposLo+1,1),plotInfo(diagix).traj(1:tposHi-tposLo+1,2), ...
            'k-','linewidth',2)
    else
        % only show start marker in individual diagram mode
        plot(plotInfo(diagix).traj(1,1),plotInfo(diagix).traj(1,2),'go','markersize',10)
        plot(plotInfo(diagix).traj(:,1),plotInfo(diagix).traj(:,2),'k-','linewidth',2)
    end
	plot(plotInfo(diagix).traj(tposHi-tposLo+1,1),plotInfo(diagix).traj(tposHi-tposLo+1,2), ...
        'ro','markersize',10)

    % start new movie for each split regime, and delete associated figure windows
    if makemovie
        if plotInfo(diagix).issplit
            switch plotInfo(diagix).frame
                case 1
                    % new movie
                    figure(fighandle(diagix))
                    frame = getframe(fighandle(diagix));
                    regime = plotInfo(diagix).regime;
                    baseFileName = [networkObjectDir '/Figures/RegimePP_' date_string '_reg' num2str(regime,'%.2d') '_000'];
                    movFileName = GetNextFilename(baseFileName, 'avi');
                    aviobj = avifile(movFileName,'colormap',colormap,'fps',10,'keyframe',5,'quality',100);
                    aviobj = addframe(aviobj,frame);
                case plotInfo(diagix).numFrames
                    % last frame
                    figure(fighandle(diagix))
                    frame = getframe(fighandle(diagix));
                    aviobj = addframe(aviobj,frame);
                    if verboseTog
                        disp(' ')
                        disp([sysStrRoot 'Saving movie file'])
                    end
                    aviobj = close(aviobj);
                    %movie2avi(frames,movFileName,'colormap',colormap,'fps',10,'keyframe',5,'quality',100)
                otherwise
                    figure(fighandle(diagix))
                    frame = getframe(fighandle(diagix));
                    aviobj = addframe(aviobj,frame);
            end
            delete(fighandle(diagix))
            savefig = false;
        else
            savefig = true;
        end
    else
        savefig = true;
    end
    if savefig
        % save figure
        regime = plotInfo(diagix).regime;
        baseFileName = [networkObjectDir '/Figures/RegimePP_' date_string '_reg' num2str(regime,'%.2d') '_diag' num2str(diagix,'%.3d') '_000'];
        figFileName = GetNextFilename(baseFileName, 'fig');
        saveas(fighandle(diagix),figFileName)
    end
end

disp(' ')
if makemovie
	disp([sysStrRoot 'Split regimes have no associated figures. Remaining figures were automatically saved'])
end
disp([sysStrRoot 'Auto-saved figures in `Figures` directory with names RegimePP_<TODAY`S DATE>_regX_diagY_Z'])
disp('         where the 3-digit Z increments to avoid overwriting old files')
disp('         (because old figures are not overwritten -- that`s your job to keep tidy!)')    
disp(' ')

% if we got this far, everything ok
returnResult = plotInfo;

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%         FUNCTIONS COPIED FROM STEPNET.m           %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function qsfp = GetQsfpVal(thisDE,varDataLine,onlyActives)
% only accepts "compiled" thisDE
global DE_NAMEi DE_GAMMA1i DE_GAMMA2i DE_ACTSWi DE_ISTAUFFILEi DE_TAURECIPi DE_GAMVARNAMEi ...
    DE_GAMVARPOWi DE_ISTGTFFILEi DE_TARGETi DE_INTVARNAMEi DE_INTVARPOWi LARGEBOUND
% V0eqn in terms of conductances and currents (for example of voltage equations):
%  curly braces imply summation
%  define {all currents} = {conductances*revpots} + {direct currents}
%  so V0 = sum {all currents} / sum {conductances}
if nargin == 2
    onlyActives = false;
end

sumG1 = SumGamma1(thisDE{DE_GAMMA1i},varDataLine,onlyActives);
if sumG1(1) ~= 0
    qsfp = ( sumG1(2) + SumGamma2(thisDE{DE_GAMMA2i},varDataLine,onlyActives) ) / sumG1(1);
else
    qsfp = LARGEBOUND;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function sum = SumGamma1(Gamma1,varDataLine,onlyActives)
% only accepts "compiled" Gamma1
% gam1term order: actSw, filefuncflag, taurecip, var, power, filefuncflag, target [, intvar, power]
global DE_NAMEi DE_GAMMA1i DE_GAMMA2i DE_ACTSWi DE_ISTAUFFILEi DE_TAURECIPi DE_GAMVARNAMEi ...
    DE_GAMVARPOWi DE_ISTGTFFILEi DE_TARGETi DE_INTVARNAMEi DE_INTVARPOWi

if nargin == 2
    onlyActives = false;
end

sum1 = 0;
sum2 = 0;
for g1t=1:length(Gamma1)
    g1term = Gamma1{g1t};
    varix = g1term{DE_GAMVARNAMEi};
    if varix == 0
        varix = g1term{DE_TARGETi};
        varVal = 1;
    else
        varVal = varDataLine(varix);
    end
    if ~onlyActives || ( onlyActives && g1term{DE_ACTSWi} )
        tau_recipVal = GetTauRVal(g1term{DE_ISTAUFFILEi},g1term{DE_TAURECIPi},...
            varVal,g1term{DE_GAMVARPOWi});
        targetVal = GetTargVal(g1term{DE_ISTGTFFILEi},g1term{DE_TARGETi},varVal,varDataLine);
        if g1term{DE_INTVARNAMEi} > 0
            sum1 = sum1 + tau_recipVal * varDataLine(g1term{DE_INTVARNAMEi})^g1term{DE_INTVARPOWi};
            sum2 = sum2 + tau_recipVal * varDataLine(g1term{DE_INTVARNAMEi})^g1term{DE_INTVARPOWi} * targetVal;
        else
            sum1 = sum1 + tau_recipVal;
            sum2 = sum2 + tau_recipVal * targetVal;
        end
    end
end
sum = [sum1, sum2];
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function sum = SumGamma2(Gamma2,varDataLine,onlyActives)
% only accepts "compiled" Gamma2
% gam2term order: actSw, filefuncflag, taurecip, var, power
global DE_NAMEi DE_GAMMA1i DE_GAMMA2i DE_ACTSWi DE_ISTAUFFILEi DE_TAURECIPi DE_GAMVARNAMEi ...
    DE_GAMVARPOWi DE_ISTGTFFILEi DE_TARGETi DE_INTVARNAMEi DE_INTVARPOWi

if nargin == 2
    onlyActives = false;
end

sum = 0;
for g2t=1:length(Gamma2)
    g2term = Gamma2{g2t};
    if ~onlyActives || ( onlyActives & g2term{DE_ACTSWi} )
        tau_recipVal = GetTauRVal(g2term{DE_ISTAUFFILEi},g2term{DE_TAURECIPi},...
            varDataLine(g2term{DE_GAMVARNAMEi}),g2term{DE_GAMVARPOWi});
        sum = sum + tau_recipVal;
    end
end
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function result = GetTargVal(sourceType,source,varVal,vLine)
switch sourceType
    case 0
		result = source;
    case 1
        result = eval( [ source '(' num2str(varVal) ')' ]);
    case 2
        result = vLine(source);
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function result = GetTauRVal(sourceType,parOrFunc,varVal,power)
if sourceType == 1
    result = eval( [ parOrFunc '(' num2str(varVal) ')' ]);
else
    result = parOrFunc * varVal^power;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% simplified version of this function from StepNet.m
function DEqns = SetActivesSwitch(DEqns, focusSet, deqnIxMap, numExt, thisAbsActs)
% do this only on COMPILED Gamma sets, for storing in either compiled or uncompiled equations
global DE_NAMEi DE_GAMMA1i DE_GAMMA2i DE_ACTSWi DE_ISTAUFFILEi DE_TAURECIPi DE_GAMVARNAMEi ...
    DE_GAMVARPOWi DE_ISTGTFFILEi DE_TARGETi DE_INTVARNAMEi DE_INTVARPOWi

for thisAbsIx = focusSet
    thisDEix = deqnIxMap(thisAbsIx);
    if thisDEix > 0 && thisAbsIx <= numExt % internal variables (even those with DE's) don't have candidate actives
        thisGamma1 = DEqns{thisDEix}{DE_GAMMA1i};
        thisGamma2 = DEqns{thisDEix}{DE_GAMMA2i};
        lenG1 = length(thisGamma1);
        lenG2 = length(thisGamma2);
        for g1t=1:lenG1
            DEqns{thisDEix}{DE_GAMMA1i}{g1t}{DE_ACTSWi} = ismember( thisGamma1{g1t}{DE_GAMVARNAMEi}, thisAbsActs );
        end
        for g2t=1:lenG2
            DEqns{thisDEix}{DE_GAMMA2i}{g2t}{DE_ACTSWi} = ismember( thisGamma2{g2t}{DE_GAMVARNAMEi}, thisAbsActs );
        end
    end % else nothing to do
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Version of that in AttEst_*.m -- utilizing fact that fast and slow vars are
% already known
function result = PrepareFocVars(fastFocVars, slowFocVars, varDataLine, varDataLine_thalf, Vix, Vval, ...
    focusSet, DEqns, allDEixMap)
% This function updates varDataLine with estimates for all fast primaryVar-dependent variables (relative to primary var)

varDataLine_temp = varDataLine;
varDataLine_temp(Vix) = Vval; % the primary focused variable
% initial values
varDataLine_est_full = varDataLine_temp;
varDataLine_est_reduced = varDataLine_temp;
% go through this once first, to freeze slow variables in reduced system
% approximation
if ~isempty(slowFocVars)
	for focVarIx=focusSet
        % don't include "self" entry, ie. focused var Vix
        if focVarIx ~= Vix
            focDEix = allDEixMap(focVarIx);
            if ismember(focVarIx,slowFocVars)
                % for slow focused vars, use constant value from half way
                % through regime for reduced system
                varDataLine_est_reduced(focVarIx) = varDataLine_thalf(focVarIx);
			end
        end
	end
end
if ~isempty(fastFocVars)
	% now approximate the fast variables
	varDataLine_est_reduced_temp = varDataLine_est_reduced;
	for focVarIx=focusSet
        % don't include "self" entry, ie. focused var Vix
        if focVarIx ~= Vix
            focDEix  = allDEixMap(focVarIx);
            if ismember(focVarIx,fastFocVars)
                % for fast focused vars, use their asymptotic value
                varDataLine_est_reduced(focVarIx) = GetQsfpVal(DEqns{focDEix},varDataLine_est_reduced_temp,true);
                varDataLine_est_full(focVarIx) = GetQsfpVal(DEqns{focDEix},varDataLine_temp,false);
            end
        end
	end
end
result.full = varDataLine_est_full;
result.reduced = varDataLine_est_reduced;
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function result = DialogBox(promptTitle,defVal,verboseTog)
% This function assumes that negative numbers will not be legal
%  values for any parameters to be got from user input!
if nargin ~=3
    disp('DialogBox:  Wrong number of arguments passed - expected 3')
    return
end
result = '';
if verboseTog
    comLineStr1 = 'Enter ';
    comLineStr2 = ' in dialog box';
    disp([ comLineStr1 lower(promptTitle) comLineStr2 ]);
end

prompt  = [ promptTitle '?' ];
ptitle  = 'Parameter entry:';
lines   = 1;

if isempty(defVal)
    def = {''};
else
    if isnumeric(defVal)
        def = {num2str(defVal)};
	elseif ischar(defVal)
        def = {defVal};
	else
        disp('DialogBox:  Internal error. `defVal` must be a numeric or string')
        return
	end
end

answer  = inputdlg(prompt,ptitle,lines,def,'off');
if ~isempty(answer)
    if ~isempty(answer{1})
        result = answer{1};
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function filename_write = GetNextFilename(filename,extStr)
% (c) Robert Clewley, BU, 2004
% Can only produce filenames for up to 999 different filenames with this extension in the local directory.
% extStr is the extension used for the filename. While this should NOT contain the '.' character,
%   this function will remove it if present in the first position of the string.
findDot = strfind(extStr,'.');
if ~isempty(findDot)
    if length(findDot) > 1
        disp(' GetNextFilename error: Too many `.` sub-extensions in extStr. These are not supported in this function')
        filename_write = '';
        return
    else
        if findDot ~=1
            disp(' GetNextFilename error: `.` was found inside extStr argument.')
            filename_write = '';
            return
        else % the dot is the first character
            extStr = extStr(2:length(extStr));
        end
    end
end

filename_write = filename; % initial value
filename_extstrip = EndStripStr(filename,'.');
trailingStr = EndRemainStr(filename_extstrip,'_',true);
foundFreeNum = false; % initial value
numCounter = 1;
while ~foundFreeNum & numCounter <= 999
	if ~isempty(trailingStr) && isNum(trailingStr)
        newNumStr = num2str(str2num(trailingStr)+numCounter,'%3.3i');
	else
        newNumStr = num2str(numCounter,'%3.3i');
	end
    filename_write = [EndStripStr(filename_extstrip,'_',true) '_' newNumStr '.' extStr];
    if ~exist(filename_write,'file')
        foundFreeNum = true;
    else
        numCounter = numCounter + 1;
    end
end

if numCounter > 999
    beep
    disp([' GetNextFilename error. Too many `.' extStr '` files (more than 999)! Cannot continue'])
    filename_write = '';
    return
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function stripStr = EndStripStr(inputStr,findChar,lastDelim,lastScope)
% strip trailing extension in a filename, using findChar to mark the beginning of the extension
if nargin==2
    lastDelim = false;
    lastScope = 5;
end
if nargin==3
    lastScope = 5;
end
if lastScope > length(inputStr)-1
    lastScope = length(inputStr)-1;
end
ix = findstr(findChar,inputStr);
if ~isempty(ix)
    if lastDelim
        if max(ix) > length(inputStr)-lastScope % filter out all except the 'lastScope' number of possibilities
                                                % (e.g. for filename extensions of known size)
            stripStr = inputStr(1:max(ix)-1); % keep all other occurrences of findChar except the final one
        else
            stripStr = inputStr;
        end
    else
        if min(ix)>1
            stripStr = inputStr(1:min(ix)-1); % keep only up to the first occurrence of findChar
        else
            stripStr = inputStr;
        end
    end
else
    stripStr = inputStr;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function stripStr = EndRemainStr(inputStr,findChar,lastDelim)
% this function returns the remainder (e.g. the extension) of a string, without the delimiter findChar 
if nargin==2
    lastDelim = false;
end
ix = findstr(findChar,inputStr);
if ~isempty(ix)
    if lastDelim
        if max(ix)>1
            stripStr = inputStr(max(ix)+1:length(inputStr));
        else
            stripStr = [];
        end
    else
        if min(ix)>=1
            stripStr = inputStr(min(ix)+1:length(inputStr));
        else
            stripStr = [];
        end
    end
else
    stripStr = [];
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function result = isNum(data)
% parameter `data` is a string
result = false;
pointCount = 0;
num = [char(48:57),'.','-','e'];
for i=1:length(data)
    if data(i) == '.'
        pointCount = pointCount + 1;
        if pointCount > 1
            return % false
        end
    end
    if ~ismember(data(i),num)
        return % false
    end
end

result = true; % can only get here if all chars were numeric
return
