function [ varargout ] = RAPTOR_plot_GUI_2012b(varargin)
%RAPTOR_PLOT_GUI
%Interactive GUI for plotting RAPTOR outputs
%
% Simple usage:
%         RAPTOR_plot_GUI(outs); % default
%         RAPTOR_plot_GUI(outs,'overview',[0 1]); % overview plot in time window
%
% Advanced usage:
% Get default parameters:
%   def = RAPTOR_plot_GUI;
%
% Get a pre-defined parameter structure
%   pp  = RAPTOR_plot_GUI(myconfig)
% where myconfig can be 'default','simple','overview'
%
% Generate the interactive figure with a customized parameter structure pp
%    hf = RAPTOR_plot_GUI(outs,pp)
% where
%    outs: a cell array of RAPTOR_out structures
%    pp  : a RAPTOR_plot_GUI parameter structure
%    hf: output figure handle
% for examples, see RAPTOR_tutorial_2 in the demos folder
%
% F. Felici TU/e 2013-2014


defplots = struct(...
    'ydata','out.Ip(end,:)',...
    'xdata','out.time',...
    'ylabel',[],...
    'xlabel',[],...
    'ylim',[],...  % optionaly axis limits
    'xlim',[],...  % optionaly axis limits
    'yscale',1,... % y scaling factor
    'axesnr',1,... % number of axis on which to plot this
    'sty','-',...
    'col',[]); % linestyle

defaxes = struct(...
    'position',[0.1 0.1 0.8 0.8]);

defstyle = struct('col',cell({'b','r',[0 0.4 0],'m',[.4 .4 .4],[0.4 0 0.4],[1 0.5 0]}));

def.plots = {defplots};
def.style = defstyle;
def.axes  = defaxes;

if nargin==0
    varargout{1} = def;
    return
elseif nargin==1
    if isstruct(varargin{1}) || iscell(varargin{1})
        RAPTOR_plot_GUI_2012b(varargin{1},'overview'); % call with default
    elseif ischar(varargin{1})
        varargout{1} = preset_config(varargin{1});
    else
        error('invalid input');
    end
    return
elseif nargin>=2
    outs = varargin{1};
    if ischar(varargin{2})
        pp = preset_config(varargin{2});
    elseif isstruct(varargin{2})
        pp = varargin{2};
    else
        error('invalid second parameter')
    end
    if nargin >= 3 && ~isempty(varargin{3})
        twin = varargin{3};
    else
        twin = []; % default, no time window
    end
    if nargin >= 4 && ~isempty(varargin{4})
        tnow = varargin{4};
    else
        tnow = []; % default, no set time
    end
end

%% checks
pp.plots = assign_defaults(pp.plots,defplots,twin);

for ip=1:numel(pp.plots)
    assert(pp.plots{ip}.axesnr<=numel(pp.axes),'axesnr exceeds number of specified axes')
end
assert(numel(outs)<=numel(pp.style),'too many out for given colors')

if numel(outs) == 1 && ~iscell(outs)
    outs = {outs}; % put single one in cell
end

%% create main plotting figure
[hf,hp,hax] = gui_init(pp,outs);
hui = ui_init(hp,hf,pp,outs);

if ~isempty(twin)
    set(hui.ht0,'string',spf(twin(1)));
    set(hui.htf,'string',spf(twin(2)));
    cb = get(hui.ht0,'callback');
    hgfeval(cb,hui.ht0); %undocumented matlab: invoke callbacks
    hgfeval(cb,hui.htf); %undocumented matlab: invoke callbacks
end
if ~isempty(tnow)
    % set current time if any
    set(hui.htime,'string',spf(tnow));
    % run callback
    cb = get(hui.htime,'callback');
    hgfeval(cb,hui.htime); %undocumented matlab: invoke callbacks
end

set(hf,'paperpositionmode','auto');

varargout{1} = hf;
varargout{2} = hax;

end

function [hf,hp,hax] = gui_init(pp,outs)

%% create and define main figure
figname = 'RAPTOR output';
hf = findall(0,'name',figname);
if numel(hf)>1; delete(hf); end
if isempty(hf) || ~ishandle(hf)
    % relative screen size
    sz = get(0,'ScreenSize'); pos = [sz(3)*0.05 sz(4)*0.2 sz(3)*0.9 sz(4)*0.7];
    % create new figure
    hf = figure('Name',figname,'NumberTitle','off','menubar','none',...
        'handlevisibility','callback','position',pos,'CloseRequestFcn',@my_closereq);
    set(hf,'units','normalized')
else % clear old figure and restart
    clf(hf); figure(hf);
end
%%
%% create axes
nax = numel(pp.axes); hax=zeros(nax,1);
for iax = 1:nax
    hax(iax) = axes('parent',hf,'position',pp.axes(iax).position);
end
hp = createplots(hax,pp,outs);

refreshdata(hf);
set(hf,'WindowButtonMotionFcn',{@wbmf,hax,hp,pp,outs});

%% nested sub-functions
    function hp = createplots(hax,pp,outs)
        %% plots
        timeplots = []; profileplots = []; % init
        
        np = numel(pp.plots);
        for iout = 1:numel(outs)
            out = outs{iout};
            ps = pp.style(iout);
            for ip=1:np
                px = pp.plots{ip};   % this plot description structure
                ax = hax(px.axesnr); % this axis handle
                
                % x data strings

                nix = numel(strfind(px.xdata,'%i'));
                if nix>0
                    xdatastr = sprintf(px.xdata,1*ones(nix,1));
                else
                    xdatastr = px.xdata;
                end
                
                niy = numel(strfind(px.ydata,'%i'));
                if niy>0
                    % this signals the profile changes in time
                    ydatastr = sprintf(px.ydata,1*ones(niy,1)); % plot first profile
                else
                    ydatastr = px.ydata;
                end
                
                % optional color override
                if ~isempty(px.col)
                    col = px.col;
                else
                    col = ps.col; % default style
                end
                
   
                
                % plot
                try
                    hpx = plot(ax,0,NaN);
                    set(hpx,'XData',eval(xdatastr),...
                        'YData',px.yscale*eval(ydatastr),...
                        'color',col,...
                        'linestyle',px.sty);
                    
                    % decide what kind of plot this is
                    if ~isempty(strfind(px.xdata,'time'))
                        timeplots = [timeplots;ax];
                    elseif ~isempty(strfind(px.xdata,false))
                        profileplots = [profileplots;ax];
                    end
                catch err
                    warning('could not plot \n Xdata:%s\n Ydata:%s\n',xdatastr,ydatastr)
                    hpx = NaN; % set an invalid handle
                end
                
                hold(ax,'on');
                
                if iout==1 % do labels only on first output
                    % x label
                    labelff(ax,px.xlabel,'x');
                    
                    % y label
                    hly = get(ax,'ylabel');
                    ystr = get(hly,'string');
                    if isempty(ystr) % generate label
                        ystr = sprintf('%s(%s)',px.ylabel,px.sty);
                    else % append new plot
                        ystr = sprintf('%s,%s(%s)',ystr,px.ylabel,px.sty);
                    end
                    labelff(ax,ystr,'y');
                end
                
                % axis limits
                if ~isempty(px.ylim)
                    set(ax,'ylim',px.ylim);
                end
                if ~isempty(px.xlim)
                    set(ax,'xlim',px.xlim);
                end
                
                % store handles
                hp(ip,iout) = hpx;
            end
        end
        if ~isempty(profileplots)
            linkaxes(profileplots,'x')
        end
        if ~isempty(timeplots)
            linkaxes(timeplots,'x')
        end
    end

    function wbmf(src,evnt,hax,hp,pp,outs)
        if strcmp(get(src,'SelectionType'),'normal')
            cp = get(src,'CurrentPoint');
            % find which axis is under pointer
            for iax=1:numel(hax);
                haxpos = get(hax(iax),'position');
                haxbox = [haxpos(1)+[0;haxpos(3)],...
                    haxpos(2)+[0;haxpos(4)]];
                
                if all(cp>=haxbox(1,:) & cp<=haxbox(2,:))
                    % find all the plots that are on this axis
                    hpxx=findobj(hp(~isnan(hp)),'parent',hax(iax));
                    if ~isempty(hpxx) % do nothing if empty
                        pxx = pp.plots(any(hp==hpxx(1),2)); % get correponding plot info for first plot
                        
                        set(src,'pointer','crosshair');
                        ax = hax(iax);
                        pt = get(ax,'Currentpoint');
                        
                        % plot values closest to tip
                        % find closet curve
                        xdata = get(hpxx,'xdata'); ydata = get(hpxx,'ydata');
                        if numel(hpxx)>1
                            % have to choose closest from multiple figures
                            for ip=1:numel(hpxx);
                                % find closest point in x
                                [~,ix(ip)] = min((xdata{ip}-pt(1)).^2);
                                yyy(ip) = abs(ydata{ip}(ix(ip))-pt(3));
                            end
                            
                            [~,ipp] = min(yyy); % find closest y
                            ipt = ix(ipp); % index of closest point
                            xx=xdata{ipp}(ipt); yy=ydata{ipp}(ipt);
                        else
                            % take closest xdata directly
                            ipp=1;
                            [~,ipt] = min((xdata-pt(1)).^2);
                            xx=xdata(ipt); yy=ydata(ipt);
                        end
                        
                        % add text at appropriate closest x,y location
                        mylabel = pp.plots{any(hp==hpxx(ipp),2)}.ylabel;
                        col = get(hpxx(ipp),'color');
                        str = sprintf('%s:%3.3g,%3.3g',mylabel,xx,yy);
                        
                        ht = findobj(src,'tag','tiptext');
                        if ~isempty(ht)
                            delete(ht); % remove if pre-existing
                        end
                        % place text string
                        text(xx,yy,str,'tag','tiptext','parent',ax,'color',col);
                        
                        if ~isempty(strfind(pxx{1}.xdata,'time'))
                            % if the pointer is inside the box axis AND
                            % the axis has 'time' on the x axis
                            
                            %fprintf('[%g,%g]\n',[pt(1),pt(3)]);
                            t_plot = pt(1); % time to plot
                            
                            % call function to change profile plot times if any
                            uidata.hp = hp; uidata.pp=pp; uidata.outs=outs; uidata.t = t_plot;
                            updateprofileplots(uidata);
                        else
                            % rho plot: do nothing
                        end
                    end
                end
            end
            
        end
    end
end

function hui = ui_init(hp,hf,pp,outs)

%%
ctrname = 'RAPTOR GUI control panel';
hc = findall(0,'name',ctrname);
if numel(hc)>1; delete(hc); end
if isempty(hc) || ~ishandle(hc)
    hc = figure('Name',ctrname,'NumberTitle','off','menubar','none','resize','off',...
        'position',[50 50 400 100],'handlevisibility','callback','CloseRequestFcn',@my_closereq);
else
    clf(hc); figure(hc); % re-create
end
%%

tmin = outs{1}.time(1);
tmax = outs{1}.time(end);
dt = median(diff(outs{1}.time));

% cosmetics
uicontrol('parent',hc,'style','text','string','spd','units','points','position',[170,40 30 20])
uicontrol('parent',hc,'style','text','string','t0' ,'units','points','position',[195,40 20 20])
uicontrol('parent',hc,'style','text','string','t'  ,'units','points','position',[230,40 20 20])
uicontrol('parent',hc,'style','text','string','tf' ,'units','points','position',[265,40 20 20])

% slider
hui.hslider = uicontrol('parent',hc,'style','slider','units','points','position',[10 10 280 20],...
    'min',tmin,'max',tmax,'value',tmin,...
    'SliderStep',[dt,0.05],'string','hello',...
    'callback',@sliderfcn);

% time indicator
hui.htime   = uicontrol('parent',hc,'style','edit','units','points','position',[220 30 40 20],...
    'string',spf(tmin),'backgroundcolor','white','tooltipstring','current time',...
    'callback',@t0tffun);

% spd
hui.hspd    = uicontrol('parent',hc,'style','edit'  ,'units','points','position',[170 30 20 20],...
    'string','1','backgroundcolor','white','tooltipstring','playback speed',...
    'callback',@setspdfcn);

% time start/stop indicator
hui.ht0      = uicontrol('parent',hc,'style','edit'  ,'units','points','position',[190 30 30 20],...
    'string',spf(tmin),'backgroundcolor','white','tooltipstring','start time',...
    'callback',@t0tffun);

hui.htf     = uicontrol('parent',hc,'style','edit'  ,'units','points','position',[260 30 30 20],...
    'string',spf(tmax),'backgroundcolor','white','tooltipstring','stop time',...
    'callback',@t0tffun);

% playback buttons
hui.hrew    = uicontrol('parent',hc,'style','pushbutton','units','points','position',[10 30 20 20],...
    'string','<<','tooltipString','Rewind',...
    'callback',@rewfcn);
hui.hrec     = uicontrol('parent',hc,'style','togglebutton','units','points','position',[70 30 20 20],...
    'string','o','tooltipString','Record movie. Toggle, then press play to start',...
    'ForegroundColor',[1 0 0],'callback','');
hui.hplay   = uicontrol('parent',hc,'style','togglebutton','units','points','position',[50 30 20 20],...
    'string','>','tooltipstring','Play','ForegroundColor',[0 0.4 0],...
    'callback',@playfcn);
hui.hstop   = uicontrol('parent',hc,'style','pushbutton','units','points','position',[30 30 20 20],...
    'string','s','tooltipstring','Stop',...
    'callback',@playfcn,'enable','off');

hui.hfaster = uicontrol('parent',hc,'style','pushbutton','units','points','position',[110 30 20 20],...
    'string','x2','tooltipstring','faster',...
    'callback',@setspdfcn);
hui.hslower = uicontrol('parent',hc,'style','pushbutton','units','points','position',[90 30 20 20],...
    'string','/2','tooltipstring','slower',...
    'callback',@setspdfcn);
hui.hloop   = uicontrol('parent',hc,'style','togglebutton','units','points','position',[130 30 40 20],...
    'string','loop','tooltipstring','loop','value',1,... % on by default
    'callback','');
%%
% plot buttons [10 30 20 20]
hui.hzoomx = uicontrol('parent',hc,'style','togglebutton','units','points','string','zoom x','position',[10  50 50 20],...
    'callback',@zoomfcn);
hui.hzoomy = uicontrol('parent',hc,'style','togglebutton','units','points','string','zoom y','position',[60  50 50 20],...
    'callback',@zoomfcn);
hui.hpan   = uicontrol('parent',hc,'style','togglebutton','units','points','string','pan'   ,'position',[110 50 30 20],...
    'callback',@zoomfcn);
%%
uidata.t  = get(hui.hslider,'value');
uidata.spd = str2double(get(hui.hspd,'string'));
uidata.t0  = str2double(get(hui.ht0,'string'));
uidata.tf  = str2double(get(hui.htf,'string'));
uidata.tmin  = tmin;
uidata.tmax  = tmax;
uidata.outs = outs;
uidata.hp = hp;
uidata.pp = pp;

guidata(hc,uidata);

% ----- nested functions have access to uidata------

    function zoomfcn(src,~)
        switch src
            case hui.hzoomx
                if get(src,'value')
                    zoom(hf,'xon');
                    set(hui.hzoomy,'value',0); % toggle other button
                    set(hui.hpan,'value',0); % toggle other button
                else
                    zoom(hf,'off');
                end
            case hui.hzoomy
                if get(src,'value')
                    zoom(hf,'yon');
                    set(hui.hzoomx,'value',0); % toggle other button
                    set(hui.hpan,'value',0); % toggle other button
                else
                    zoom(hf,'off');
                end
            case hui.hpan
                if get(src,'value')
                    pan(hf,'on');
                    set(hui.hzoomx,'value',0); % toggle other button
                    set(hui.hzoomy,'value',0); % toggle other button
                else
                    pan(hf,'off');
                end
        end
    end

    function sliderfcn(src,~)
        t=get(src,'value'); % get time on slider
        set(hui.htime,'string',spf(t));
        
        uidata=guidata(src); uidata.t=t;
        if t>uidata.tf
            uidata.tf=t;
            set(hui.htf,'string',spf(t));
        elseif t<uidata.t0
            uidata.t0=t;
            set(hui.ht0,'string',spf(t));
        end
        
        setaxislims
        updateprofileplots(uidata);
        guidata(src,uidata);
        
    end

    function rewfcn(src,~)
        uidata=guidata(src);
        switch src
            case hui.hrew % rewind
                t = uidata.t0;
            case hui.hff % forward
                t = uidata.tf;
        end
        set(hui.hslider,'value',t);
        set(hui.htime,'string',spf(t));
        
        uidata.t=t;
        updateprofileplots(uidata);
        guidata(hc,uidata);
    end

    function setspdfcn(src,~)
        currentspd = str2double(get(hui.hspd,'string'));
        
        switch src
            case hui.hslower
                set(hui.hspd,'string',num2str(currentspd/2));
            case hui.hfaster
                set(hui.hspd,'string',num2str(currentspd*2));
            case hui.spd % called from text box directly
                
            otherwise
                error('unknown mode')
        end
        
    end

    function playfcn(src,~)
        
        switch src
            case hui.hplay
                tic;
                if get(src,'value')==1 % play
                    % check that we are not at tmax already
                    tf = str2double(get(hui.htf,'string'));
                    t0 = str2double(get(hui.ht0,'string'));

                    t_plot = get(hui.hslider,'value'); % start value for plot
                    
                    if t_plot<=tf
                        set(hui.hstop,'enable','on');
                        set(src,'string','||','tooltipstring','Pause');
                        
                        dorec = get(hui.hrec,'value');
                        set(hui.hrec,'enable','off');
                        
                        if dorec % preparations for recording phase
                            t_plot = t0; % start time
                            set(hui.hslider,'value',t0); % go to t0 to start
                            set(hui.hloop,'value',0); % disable looping
                            % init recording
                            fprintf('recording movie for t=[%3.3f,%3.3f]\n',t0,tf); 
                            [writerObj,moviefile] = init_recording;
                            if isempty(writerObj)
                                fprintf('failed to make writerObj. Recording aborted. \n');
                                set(hui.hrec,'value',0); % stop recording
                                set(src,'value',0,'string','>','tooltipstring','Play'); % stop playing
                                dorec = 0;
                            end
                            set(hf,'Renderer','zbuffer');
                        end
                        
                    end
                    
                    while get(src,'value')==1
                        tic; % timer
                        % reached limit
                        tf = str2double(get(hui.htf,'string'));
                        if t_plot>tf;
                            doloop = get(hui.hloop,'value');
                            if doloop % loop: start again
                                t_plot = str2double(get(hui.ht0,'string'));
                            else
                                % stop here
                                
                                % set last point as current time
                                uidata=guidata(src);
                                uidata.t=tf;
                                guidata(hc,uidata); % save
                                set(hui.htime,'string',num2str(tf));

                                set(src,'value',0);
                                playfcn(src,[]);
                                break
                            end
                        end
                        set(hui.hslider,'value',t_plot);
                        set(hui.htime,'string',num2str(t_plot));
                        
                        % set new time
                        uidata=guidata(src); uidata.t=t_plot;
                        % update with new data
                        updateprofileplots(uidata);
                        guidata(hc,uidata); % save
                        
                        spd = str2double(get(hui.hspd,'string'));
                        if ~dorec
                           % next plot time based on elapsed time
                            t_plot = t_plot + spd*(toc);
                        else % if recording
                            % save frame
                            perc = (t_plot-t0)/(tf-t0)*100;
                            fprintf('%02.f%%, t=%3.3f\n',perc,t_plot); % progress
                            frame = getframe(hf);
                            writeVideo(writerObj,frame);
                            % advance time at 24fps * speed
                            fps = 24;
                            t_plot = t_plot + spd/fps;
                        end
                    end
                    % re-enable stuff
                    
                    if dorec % close rec object
                       close(writerObj);
                       fprintf('Finished writing movie to: %s\n',moviefile);
                       set(hui.hrec,'value',0)
                    end
                    set(hui.hrec,'enable','on')
                else
                    set(src,'string','>','tooltipstring','Play');
                end
            case hui.hstop
                set(hui.hplay,'value',0,'string','>'); % stop playback
                
                % reset to first time
                t=uidata.t0;
                set(hui.hslider,'value',t);
                set(hui.htime,'string',spf(t));
                uidata.t=t;
                updateprofileplots(uidata);
                guidata(hc,uidata);
                set(hui.hstop,'enable','off');
        end
        
    end

    function t0tffun(src,~) % callback when setting time limits
        uidata=guidata(src);
        switch src
            case hui.ht0
                t0set = str2double(get(hui.ht0,'string'));
                if ~isnumeric(t0set) || t0set<uidata.tmin || t0set>uidata.tf
                    % if invalid reset old
                    warning('invalid setting for t0')
                    set(hui.ht0,'string',spf(uidata.t0));
                else % set new limit
                    uidata.t0 = t0set;
                    if uidata.t<uidata.t0
                        uidata.t=uidata.t0;
                    end
                end
                
            case hui.htf
                tfset = str2double(get(hui.htf,'string'));
                if ~isnumeric(tfset) || tfset>uidata.tmax || tfset<uidata.t0
                    % if invalid reset old
                    warning('invalid setting for tf')
                    set(hui.htf,'string',spf(uidata.tf));
                else % else set new limit
                    uidata.tf = tfset;
                    if uidata.t>uidata.tf
                        uidata.t = uidata.tf;
                    end
                end
                
            case hui.htime
                t = str2double(get(hui.htime,'string'));
                if ~isnumeric(t) || t>uidata.tmax || t<uidata.tmin
                    warning('invalid setting for t')
                    set(hui.htime,'string',spf(uidata.t));
                else % set new time
                    uidata.t = t;
                    if t>uidata.tf
                        uidata.tf=t;
                        set(hui.htf,'string',spf(t));
                    elseif t<uidata.t0
                        uidata.t0=t;
                        set(hui.ht0,'string',spf(t));
                    end
                end
        end
        
        set(hui.hslider,'value',uidata.t);
        set(hui.htime,'string',spf(uidata.t));
        updateprofileplots(uidata);
        guidata(hc,uidata); % save
        
        setaxislims; % set new axis limits
    end

    function setaxislims
        for iout = 1:numel(outs)
            for ip=1:numel(pp.plots)
                hpx = hp(ip,iout);
                px = pp.plots{ip};
                
                if isnan(hpx) || ~isempty(strfind(px.xdata,false)) % if NaN or rho on x axis
                    % do nothing
                elseif ~isempty(strfind(px.xdata,'time')) % if time on x axis
                    % change axis limits
                    ax = get(hpx,'parent'); % get parent axis
                    
                    set(ax,'xlim',[uidata.t0,uidata.tf])
                end
            end
        end
    end
end

function updateprofileplots(uidata)

hp=uidata.hp; pp=uidata.pp; outs=uidata.outs; t_plot=uidata.t;

for iout = 1:numel(outs)
    out = outs{iout};
    
    % find time index to plot
    it = find(out.time>=t_plot,1,'first');
    if isempty(it)
        it=numel(out.time);
    end
    
    % plot
    for ip=1:numel(pp.plots)
        hpx = hp(ip,iout);
        px = pp.plots{ip};
        
        if ishandle(hpx) % only if it is a valid handle
            % add it to string if necessary
            niy = numel(strfind(px.ydata,'%i'));
            if niy>0
                try
                    ydatastr = sprintf(px.ydata,it*ones(1,niy));
                    set(hpx,'ydata',px.yscale*eval(ydatastr));
                catch err
                    % ignore
                end
            end
            
            nix = numel(strfind(px.xdata,'%i'));
            if nix>0
                try
                    xdatastr = sprintf(px.xdata,it*ones(1,nix));
                    set(hpx,'xdata',eval(xdatastr));
                catch err
                    % ignore
                end
            end
            
            % draw vertical line if necessary
            if ~isempty(strfind(px.xdata,'time')) % if time on x axis
                % draw vertical line
                ax = get(hpx,'parent'); % get parent axis
                hv = findobj(ax,'tag','vertical line'); % find an object with this tag
                
                % x and y values for vertical line
                yy = get(ax,'ylim');
                xx = [t_plot,t_plot];
                
                % plot new object, or change old object dpeending
                if isempty(hv)
                    hv = plot(ax,xx,yy,'color',0.5*[1 1 1]);
                    set(hv,'tag','vertical line');
                else
                    set(hv,'xdata',xx,'ydata',yy,'linestyle','--','linewidth',1);
                end
            end
        end
    end
end
drawnow

end

function my_closereq(src,evnt)
% close request function that closes both windows (window and control panel) if
% either is closed

ctrname = 'RAPTOR GUI control panel';
hc = findall(0,'name',ctrname);

figname = 'RAPTOR output';
hf = findall(0,'name',figname);

delete(hf);
delete(hc);
end

function hl = labelff(varargin)
% h = labelff(string,xy)
% or
% h = labelff(axishandle,string,xy)
% puts x or y label "string" inside axes in a corner.
% xy = 'x' or 'y'

if nargin==1;
    error('insufficient inputs')
elseif nargin>1
    if ishandle(varargin{1})
        hax = varargin{1}; % axis handle
        string = varargin{2};
        xy = varargin{3};
    else
        hax = gca;
        string = varargin{1};
        xy = varargin{2};
    end
end

if isnumeric(xy)
    assert(numel(xy)==2,'xy must be x,y, or 2 numbers');
    pos = xy;
    if xy(1)>xy(2)
        hl = xlabel(hax,string);
    else
        hl = ylabel(hax,string);
    end
else
    switch xy
        case 'x'
            hl = xlabel(hax,string);
            pos = [0.99,0.01];
            set(hl,'units','normalized','position',pos,'horizontalalignment','right','verticalalignment','bot','rotation',0)
            
        case 'y'
            hl = ylabel(hax,string);
            pos = [0.01 0.97];
            set(hl,'units','normalized','position',pos,'horizontalalignment','left','verticalalignment','top','rotation',0)
        otherwise;
            error('choose x or y')
    end
end
end

function pp_mod = assign_defaults(pp,def,twin)

pp_mod = pp;
fields = fieldnames(def);

for ip=1:numel(pp)
    
    for ii=1:numel(fields)
        field = fields{ii};
        % add default to missing field
        if ~isfield(pp{ip},field) || isempty(pp{ip}.(field))
            pp_mod{ip}.(field) = def.(field);
        end
        % add time window data to time plots
        if ~isempty(twin) && ~isempty(strfind(pp{ip}.xdata,'time'))
            pp_mod{ip}.xlim = twin;
        end
    end
end

end

function positions = axespos(nrows,ncols,rowmarg,colmarg,spaces)

rowleft = rowmarg(1);
rowright = rowmarg(2);
coltop = colmarg(1);
colbot = colmarg(2);
srows = spaces(1);
scols = spaces(2);

drows = ((1-rowleft-rowright) - srows*(nrows-1))/nrows;
dcols = ((1-coltop-colbot) - scols*(ncols-1))/ncols;

for irow = 1:nrows;
    for icol = 1:ncols;
        prow = rowleft+(drows+srows)*(irow-1);
        pcol = 1-coltop-icol*dcols - scols*(icol-1);
        positions(irow,icol).position = [prow,pcol,drows,dcols];
    end
end
end

function str=spf(num)
str = sprintf('%3.3g',num);
end

function [writerObj,moviefile] = init_recording
% initialize recording

pathstr = fileparts(mfilename('fullpath'));
moviepath = fullfile(pathstr,'..','personal','movies','*.mp4');
[filename,pathname] = uiputfile(moviepath);

% if user aborts here, or other problem, return empty and abort recording
if ~ischar(filename)
    writerObj = [];
    moviefile = ''; 
    return
else
    % prepare movie object
    moviefile = fullfile(pathname,filename);
    writerObj = VideoWriter(moviefile,'Motion JPEG AVI');
    set(writerObj,'FrameRate',24);
    open(writerObj);
end

end

function pp = preset_config(myconfig)

switch myconfig
    case 'default'
        pp = RAPTOR_plot_GUI_2012b; % call with no input returns def
    case 'simple'
        pp = RAPTOR_plot_GUI_2012b; % get defaults
        
        % simple plots with some key quantities
        p{1} = struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ip(end,:)' ,'ylabel','I_p[MA]','yscale',1e-6,'axesnr',1,'sty','-');
        
        p{2} = struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.q(:,%i)' ,'ylabel','q','yscale',1,'axesnr',2,'sty','-');
        
        p{3} = struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.te(:,%i)' ,'ylabel','T_e[keV]','yscale',1e-3,'axesnr',3,'sty','-');
        
        p{4} = struct('xdata','out.time','xlabel','t',...
            'ydata','out.te(1,:)' ,'ylabel','T_{e0}[keV]','yscale',1e-3,'axesnr',4,'sty','-');
        
        p{5} = struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ne(:,%i)' ,'ylabel','n_{e19}','yscale',1e-19','axesnr',5,'sty','-');
        
        p{6} = struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ti(:,%i)' ,'ylabel','T_i[keV]','yscale',1e-3','axesnr',3,'sty','--');
        
        %         p(1).position = [0.1 0.7 0.3 0.2];
        %         p(2).position = [0.1 0.4 0.3,0.2];
        %         p(3).position = [0.1 0.1 0.3 0.2];
        %         p(4).position = [0.5 0.7 0.3 0.2];
        %         p(5).position = [0.5 0.4 0.3,0.2];
        %         p(6).position = [0.5 0.1 0.3 0.2];
        
        pp.plots = p; % reassign
        
        % axes position definitions
        pp.axes = axespos(2,3,[0.05 0.05],[0.05 0.05],[0.05,0.05]);
        
    case 'hcd' % key quantities related to heating and current drive
        pp = RAPTOR_plot_GUI_2012b; % get defaults
        
        p={}; % init
        % Ip
        p = [p,{struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ip(end,:)' ,'ylabel','I_p[MA]','yscale',1e-6,'axesnr',1,'sty','-')}];
        
        % loop voltage
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.upl(end,:)' ,'ylabel','V_{loop}[V]','yscale',1,'axesnr',1,'sty','-')];
        
        % internal inductance
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.li3(end,:)' ,'ylabel','li3','yscale',1,'axesnr',1,'sty','--')];
        
        % q profile info
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.qmin' ,'ylabel','q_{min}','yscale',1,'axesnr',2,'sty','-')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.q0' ,'ylabel','q_0','yscale',1,'axesnr',2,'sty',':')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.q95' ,'ylabel','q_{95}','yscale',1,'axesnr',2,'sty','--')];
        
        % temperature and density
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.te(1,:)' ,'ylabel','T_{e0}[keV]','yscale',1e-3,'axesnr',3,'sty','-')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.ne(1,:)' ,'ylabel','n_{e0}[10^{19}m^{-3}]','yscale',1e-19,'axesnr',3,'sty','--')];
        
        % input powers
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Pec(end,:)' ,'ylabel','\SigmaP_{ec}','yscale',1e-6,'axesnr',4,'sty','--')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Pnb(end,:)' ,'ylabel','\SigmaP_{nb}[MW]','yscale',1e-6,'axesnr',4,'sty','-.')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Poh(end,:)' ,'ylabel','\SigmaP_{oh}','yscale',1e-6,'axesnr',4,'sty','-')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Palphatot(end,:)' ,'ylabel','\SigmaP_{\alpha}[MA]','yscale',1e-6,'axesnr',4,'sty',':')];
        
        % current components
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ioh(end,:)' ,'ylabel','I_{oh}','yscale',1e-6,'axesnr',5,'sty','-')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ibs(end,:)' ,'ylabel','I_{bs}','yscale',1e-6,'axesnr',5,'sty',':')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Inb(end,:)' ,'ylabel','I_{nb}','yscale',1e-6,'axesnr',5,'sty','--')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Iec(end,:)' ,'ylabel','I_{ec}','yscale',1e-6,'axesnr',5,'sty','-.')];
        
        % Upl
        p = [p, struct(...
            'xdata','out.rho',      'xlabel','\rho',...
            'ydata','out.upl(:,%i)','ylabel','U_{pl}[V]','yscale',1,'axesnr',6,'sty','-')];
        
        % q and shear
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.q(:,%i)' ,'ylabel','q','yscale',1,'axesnr',6,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.shear(:,%i)' ,'ylabel','shear','yscale',1,'axesnr',6,'sty','--')];
        % rational surfaces
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','[1]*ones(size(out.rho))' ,'ylabel','','yscale',1,'axesnr',6,'sty',':')];
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','[2]*ones(size(out.rho))' ,'ylabel','','yscale',1,'axesnr',6,'sty',':')];
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','[3]*ones(size(out.rho))' ,'ylabel','','yscale',1,'axesnr',6,'sty',':')];
        
        
        % driven currents
        %p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %    'ydata','out.jaux(:,%i)' ,'ylabel','j_{cd}','yscale',1e-6,'axesnr',14,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.jec(:,%i)' ,'ylabel','j_{ec}','yscale',1e-6,'axesnr',8,'sty','--')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.jnb(:,%i)' ,'ylabel','j_{nb}','yscale',1e-6,'axesnr',8,'sty','-.')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.jbs(:,%i)' ,'ylabel','j_{bs}[MA/m^2]','yscale',1e-6,'axesnr',8,'sty',':')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.joh(:,%i)' ,'ylabel','j_{oh}[MA/m^2]','yscale',1e-6,'axesnr',8,'sty','--')];
      
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.jpar(:,%i)' ,'ylabel','j_{||}[MA/m^2]','yscale',1e-6,'axesnr',8,'sty','-')];
      
        
        %         % cumulative current components
        %         p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %             'ydata','out.joh(:,%i)' ,'ylabel','j_{\Omega}','yscale',1e-6,'axesnr',15,'sty','-')];
        %
        %         callstr = 'subsref(out.joh+out.jbs,struct(''type'',''()'',''subs'',{{'':'',%i}}))';
        %         p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %             'ydata', callstr ,'ylabel','+j_{bs}','yscale',1e-6,'axesnr',15,'sty','--')];
        %
        %         callstr = 'subsref(out.joh+out.jbs+out.jaux,struct(''type'',''()'',''subs'',{{'':'',%i}}))';
        %         p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %             'ydata',callstr,'ylabel','+j_{aux}[MA/m^2]','yscale',1e-6,'axesnr',15,'sty','-')];
        %
        
        % kinetic profiles
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.te(:,%i)' ,'ylabel','T_e[keV]','yscale',1e-3,'axesnr',9,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ti(:,%i)' ,'ylabel','T_i[keV]','yscale',1e-3,'axesnr',9,'sty','--')];
        
        
        % density
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ne(:,%i)' ,'ylabel','n_{e}[10^{19}m^{-3}]','yscale',1e-19,'axesnr',9,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ni(:,%i)' ,'ylabel','n_{i}[10^{19}m^{-3}]','yscale',1e-19,'axesnr',9,'sty',':')];
        
        % power density components
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.poh(:,%i)' ,'ylabel','p_{\Omega}','yscale',1e-6,'axesnr',7,'sty','-')];
        
        callstr = 'out.palphae(:,%i)';
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata', callstr ,'ylabel','p_{\alpha}','yscale',1e-6,'axesnr',7,'sty',':')];
        
        callstr = 'out.pnbe(:,%i)';
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata',callstr,'ylabel','p_{nbe}[MW/m^3]','yscale',1e-6,'axesnr',7,'sty','-.')];
        
        callstr = 'out.pec(:,%i)';
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata',callstr,'ylabel','p_{ec}[MW/m^3]','yscale',1e-6,'axesnr',7,'sty','--')];
        
        % automatically spaced positions for axes
        pp.plots = p; % reassign
        pp.axes =  axespos(3,3,[0.05 0.05],[0.05 0.05],[0.05 0.05]);
        
        %%%% OVERVIEW CASE %%%%%
    case 'overview' % key quantities related to current diffusion
        pp = RAPTOR_plot_GUI_2012b; % get defaults
        
        p = {}; % init
        % Ip
        p = [p,{struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ip(end,:)' ,'ylabel','I_p[MA]','yscale',1e-6,'axesnr',1,'sty','-')}];
        
        % Q
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.Q' ,'ylabel','Q_{fus}','yscale',1,'axesnr',1,'sty','-')];
        
        % q profile info
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.qmin' ,'ylabel','q_{min}','yscale',1,'axesnr',2,'sty','-')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.q0' ,'ylabel','q_0','yscale',1,'axesnr',2,'sty',':')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.q95' ,'ylabel','q_{95}','yscale',1,'axesnr',2,'sty','--')];
        
        % loop voltage
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.upl(end,:)' ,'ylabel','V_{loop}[V]','yscale',1,'axesnr',3,'sty','-')];
        
        % internal inductance
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.li3(end,:)' ,'ylabel','li3','yscale',1,'axesnr',3,'sty','--')];
        
        
        % temperature and density
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.te(1,:)' ,'ylabel','T_{e0}[keV]','yscale',1e-3,'axesnr',4,'sty','-')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.ne(1,:)' ,'ylabel','n_{e0}[10^{19}m^{-3}]','yscale',1e-19,'axesnr',4,'sty','--')];
        
        % Figures of merit
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.H98' ,'ylabel','H_{98y3}','yscale',1,'axesnr',5,'sty','--')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.tauE' ,'ylabel','tau_E','yscale',1,'axesnr',5,'sty','--')];
        
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.psi(end,:)-out.psi(end,1)' ,'ylabel','\psi_b','yscale',1,'axesnr',5,'sty',':')];
        
        
        % input powers
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Pec(end,:)' ,'ylabel','\SigmaP_{ec}','yscale',1e-6,'axesnr',6,'sty','--')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Pnb(end,:)' ,'ylabel','\SigmaP_{nbe}[MW]','yscale',1e-6,'axesnr',6,'sty','-.')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Poh(end,:)' ,'ylabel','\SigmaP_{oh}','yscale',1e-6,'axesnr',6,'sty','-')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Palphatot(end,:)' ,'ylabel','\SigmaP_{\alpha}[MA]','yscale',1e-6,'axesnr',6,'sty',':')];
        
        % Energies
        % total thermal energy
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.Wth' ,'ylabel','W_{th}','yscale',1e-6,'axesnr',7,'sty',':')];
        
        % electron energy
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.We' ,'ylabel','W_{e}','yscale',1e-6,'axesnr',7,'sty','-')];
        
        % ion energy
        p = [p,struct('xdata','out.time','xlabel','t',...
            'ydata','out.Wi' ,'ylabel','W_{i}','yscale',1e-6,'axesnr',7,'sty','--')];
        
        % power losses
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Prad(end,:)' ,'ylabel','P_{rad}','yscale',1e-6,'axesnr',8,'sty','-')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Pbrem(end,:)' ,'ylabel','P_{brem}','yscale',1e-6,'axesnr',8,'sty','--')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Pei(end,:)' ,'ylabel','P_{ei}[MW]','yscale',1e-6,'axesnr',8,'sty','-.')];
        
        % current components
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ioh(end,:)' ,'ylabel','I_{oh}','yscale',1e-6,'axesnr',9,'sty','-')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ibs(end,:)' ,'ylabel','I_{bs}','yscale',1e-6,'axesnr',9,'sty',':')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Inb(end,:)' ,'ylabel','I_{nb}','yscale',1e-6,'axesnr',9,'sty','--')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.Iec(end,:)' ,'ylabel','I_{ec}','yscale',1e-6,'axesnr',9,'sty','-.')];
        
        
        % beta
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.beta' ,'ylabel','\beta','yscale',1,'axesnr',10,'sty','-')];
        
        % p = [p, struct('xdata','out.time','xlabel','t',...
        %     'ydata','out.betap' ,'ylabel','\beta_p','yscale',1,'axesnr',10,'sty','-.')];
        
        p = [p, struct('xdata','out.time','xlabel','t',...
            'ydata','out.betaN' ,'ylabel','\beta_N','yscale',1,'axesnr',10,'sty','--')];
        
        % profiles
        
        % Upl
        p = [p, struct(...
            'xdata','out.rho',      'xlabel','\rho',...
            'ydata','out.upl(:,%i)','ylabel','U_{pl}[V]','yscale',1,'axesnr',11,'sty','-')];
        
        % q and shear
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.q(:,%i)' ,'ylabel','q','yscale',1,'axesnr',12,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.shear(:,%i)' ,'ylabel','shear','yscale',1,'axesnr',12,'sty','--')];
        
        % rational surfaces
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','[1]*ones(size(out.rho))' ,'ylabel','','yscale',1,'axesnr',12,'sty',':')];
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','[2]*ones(size(out.rho))' ,'ylabel','','yscale',1,'axesnr',12,'sty',':')];
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','[3]*ones(size(out.rho))' ,'ylabel','','yscale',1,'axesnr',12,'sty',':')];
        
        
        % power loss distributions
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.prad(:,%i)' ,'ylabel','p_{rad}[MW/m^3]','yscale',1e-6,'axesnr',13,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.pbrem(:,%i)' ,'ylabel','p_{brem}[MW/m^3]','yscale',1e-6,'axesnr',13,'sty','--')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.pei(:,%i)' ,'ylabel','p_{ei}[MW/m^3]','yscale',1e-6,'axesnr',13,'sty','-.')];
        
        % driven currents
        %p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %    'ydata','out.jaux(:,%i)' ,'ylabel','j_{cd}','yscale',1e-6,'axesnr',14,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.jec(:,%i)' ,'ylabel','j_{ec}','yscale',1e-6,'axesnr',14,'sty','--')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.jnb(:,%i)' ,'ylabel','j_{nb}','yscale',1e-6,'axesnr',14,'sty','-.')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.jbs(:,%i)' ,'ylabel','j_{bs}[MA/m^2]','yscale',1e-6,'axesnr',14,'sty',':')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.joh(:,%i)' ,'ylabel','j_{oh}[MA/m^2]','yscale',1e-6,'axesnr',14,'sty','-')];
        
        
        %         % cumulative current components
        %         p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %             'ydata','out.joh(:,%i)' ,'ylabel','j_{\Omega}','yscale',1e-6,'axesnr',15,'sty','-')];
        %
        %         callstr = 'subsref(out.joh+out.jbs,struct(''type'',''()'',''subs'',{{'':'',%i}}))';
        %         p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %             'ydata', callstr ,'ylabel','+j_{bs}','yscale',1e-6,'axesnr',15,'sty','--')];
        %
        %         callstr = 'subsref(out.joh+out.jbs+out.jaux,struct(''type'',''()'',''subs'',{{'':'',%i}}))';
        %         p = [p, struct('xdata','out.rho','xlabel','\rho',...
        %             'ydata',callstr,'ylabel','+j_{aux}[MA/m^2]','yscale',1e-6,'axesnr',15,'sty','-')];
        %
        
        % zeff
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ze(:,%i)' ,'ylabel','Z_{eff}','yscale',1,'axesnr',15,'sty','-')];
        
        % conductivity
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.signeo(:,%i)' ,'ylabel','\sigma_{||}[m\Omegam]','yscale',1e-6,'axesnr',15,'sty','-')];
        
        % kinetic profiles
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.te(:,%i)' ,'ylabel','T_e[keV]','yscale',1e-3,'axesnr',16,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ti(:,%i)' ,'ylabel','T_i[keV]','yscale',1e-3,'axesnr',16,'sty','--')];
        
        
        % iota profile
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.iota(:,%i)' ,'ylabel','\iota=1/q','yscale',1,'axesnr',17,'sty','-')];
        
        
        % density
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ne(:,%i)' ,'ylabel','n_{e}[10^{19}m^{-3}]','yscale',1e-19,'axesnr',18,'sty','-')];
        
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.ni(:,%i)' ,'ylabel','n_{i}[10^{19}m^{-3}]','yscale',1e-19,'axesnr',18,'sty',':')];
        
        % power density components
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.poh(:,%i)' ,'ylabel','p_{\Omega}','yscale',1e-6,'axesnr',19,'sty','-')];
        
        callstr = 'out.palphae(:,%i)';
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata', callstr ,'ylabel','p_{\alpha}','yscale',1e-6,'axesnr',19,'sty',':')];
        
        callstr = 'out.pnbe(:,%i)';
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata',callstr,'ylabel','p_{nbe}[MW/m^3]','yscale',1e-6,'axesnr',19,'sty','-.')];
        
        callstr = 'out.pec(:,%i)';
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata',callstr,'ylabel','p_{ec}[MW/m^3]','yscale',1e-6,'axesnr',19,'sty','--')];
        
        % chi
        p = [p, struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.chie(:,%i)' ,'ylabel','\chi_e[m^2/s]','yscale',1,'axesnr',20,'sty','-')];
         
        % automatically spaced positions for axes
        pp.plots = p; % reassign
        pp.axes =  axespos(5,4,[0.05 0.05],[0.05 0.05],[0.05 0.05]);
        
        
    case 'geometry'
        % plots related to geometry
        pp = RAPTOR_plot_GUI_2012b; % get defaults
        
        % simple plots with some key quantities
        p = {}; % init
        p = [p,{struct('xdata','out.time','xlabel','t',...
            'ydata','out.Ip(end,:)' ,'ylabel','I_p[MA]','yscale',1e-6,'axesnr',1,'sty','-')}];
        
        p = [p,struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.g1(:,%i)' ,'ylabel','g1','yscale',1,'axesnr',2,'sty','-')];
        
        p = [p,struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.g2(:,%i)' ,'ylabel','g2','yscale',1,'axesnr',3,'sty','-')];
        
        p = [p,struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.g3(:,%i)' ,'ylabel','g3','yscale',1,'axesnr',4,'sty','-')];
        
        p = [p,struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.Vp(:,%i)' ,'ylabel','Vprime','yscale',1','axesnr',5,'sty','-');];
                
        p = [p,struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.F(:,%i)' ,'ylabel','F=RB_\phi','yscale',1','axesnr',6,'sty','-');];
        
        p = [p,struct('xdata','out.rho','xlabel','\rho',...
            'ydata','out.g23or(:,%i)','ylabel','g_2g_3/\rho','yscale',1','axesnr',9,'sty','-');];
        
        p = [p,struct('xdata','out.time','xlabel','\rho',...
            'ydata','out.Phib','ylabel','\Phib','yscale',1','axesnr',7,'sty','-');];
        
        p = [p,struct('xdata','out.time','xlabel','\rho',...
            'ydata','out.Volume(end,:)','ylabel','Volume','yscale',1','axesnr',8,'sty','-');];
       
        p = [p,struct('xdata','out.time','xlabel','\rho',...
            'ydata','out.g23or(end,:)/100','ylabel','g2g3/r(end)','yscale',1','axesnr',8,'sty','--');];

        
        pp.plots = p; % reassign
        
        % axes position definitions
        pp.axes = axespos(3,3,[0.05 0.05],[0.05 0.05],[0.05,0.05]);
        
    case 'energy' % key quantities related to energy transport
    otherwise
        error('invalid option, use ''default'' or ''simple'' or ''overview''')
end
end

