function VMS(arg,options)
%VMS(arg,options)

    global gMISC_GLOBALS
    %     global VMSfig
    if nargin<2
        options='';
    end
    if nargin<1
%         GUI=get(0,'UserData');
%         if ~isfield(GUI,'VMS')
%             clear
%         if exist('VMSfig')
%             if ishandle(VMSfig)
%                 clear VMSfig
%             end
%         end
        arg='Initialise';
        data=Initialise(arg,options);
        VMSfig=data.handles.figMain;
        GUI=get(0,'UserData');
        GUI.VMS=data;
        set(0,'UserData',GUI);
        %         guidata(data.handles.figMain,data);
    end
    GUI=get(0,'UserData');
    data=GUI.VMS;
    %     data=guidata(VMSfig);
    currentfig=gcf;
    set(currentfig,'Pointer','watch')
    if ~strcmpi(arg,'Quit')
        switch arg
            case 'Initialise'
                % already initialised
            case 'Play'
                data=Play(data,options)
            case 'Project'
                data=Project(data,options);
            case 'Progress'
                data=Progress(data,options);
            case 'VariableDisplayed'
                data=VariableDisplayed(data,1,options);
            case 'Movie'
                data=Movie(data,options);
            case 'StartNumber'
                data=Startnumber(data,options);
            otherwise
                error('VMS: unrecognised argument')
        end
        %         guidata(data.handles.figMain,data);
        GUI.VMS=data;
        set(0,'UserData',GUI);
        set(currentfig,'Pointer','arrow')
    else
        cd(data.startdirectory);
        for i=1:length(data.current_figure)
            if data.current_figure(i)~=0 && ishandle(data.current_figure(i))
                close(data.current_figure(i));
            end
        end
        close(data.handles.figMain);
        GUI=rmfield(GUI,'VMS');
        set(0,'UserData',GUI);
    end
end

function data=Startnumber(data,options)
    if options==1 %start
        v=get(data.handles.Progress,'Value');
        n=ceil(v*max(data.meshnumbers));
        dist=(data.meshnumbers-n).^2;
        [dists,ii]=min(dist);
        k=ii;
        n=str2num(get(data.handles.StartNumber,'String'));;
        if ~isempty(n)
            k=n;
        end
        set(data.handles.StartNumber,'String',num2str(k));   
        data.StartNumber=k;
        v=k/max(data.meshnumbers);
        set(data.handles.Progress,'Value',v);
    else
        k=ceil(max(data.meshnumbers));
        n=str2num(get(data.handles.EndNumber,'String'));;
%        n=input(sprintf('Finish frame number [%d] ',k));
        if ~isempty(n)
            if n<=k
                k=n;
            end
        end
        set(data.handles.EndNumber,'String',num2str(k));        
        data.EndNumber=k;
    end
end

function data=Movie(data,options)
    currentd=pwd;
    if get(data.handles.Movie,'Value')
        datetime=['_',datestr(now,31)];
        datetime(findstr(datetime,' '))='';
        datetime(findstr(datetime,':'))='';
        datetime(findstr(datetime,'-'))='';
        for i=1:size(data.mdir,1)
%             tempname=fullfile(data.mdir{i},['movie',datetime,'.avi']);
            cd(data.mdir{i});
            data.moviename{i}=['movie',datetime,'.avi'];
            data.movieflags(i)=1;
            data.mov{i} = avifile(data.moviename{i},'fps',5,'compression','None','Quality',90);
        end
    else
        set(data.handles.Movie,'Value',0)
        for i=1:size(data.mdir,1)
            data.mov{i}=close(data.mov{i});
            data.movieflags(i)=0;
            disp(sprintf('Written movie to %s',data.moviename{i}));
            cd(data.mdir{i});
            nn=data.moviename{i};
            name=[nn(1:end-3),'fig'];
            saveas(data.current_axes,name);
        end
    end
    cd(currentd);
end

% m.plotdefaults.autoScale
% m.plotdefaults.axisRange

function data=Initialise(arg,options)
    fig = openfig(mfilename, 'reuse');
    movegui(fig,'northwest');
    handles = guihandles(fig);
    data.handles=handles;
    set(fig, 'Visible', 'on');
    data.handles.figMain = fig;
    set(data.handles.Project1,'callback',[mfilename,'(''Project'',1);']);
    set(data.handles.Project2,'callback',[mfilename,'(''Project'',2);'],'Visible','off');
    set(data.handles.Project3,'callback',[mfilename,'(''Project'',3);'],'Visible','off');
    set(data.handles.Play,'callback',[mfilename,'(''Play'',1);']);
    set(data.handles.Progress,'callback',[mfilename,'(''Progress'');']);
    set(data.handles.End,'callback',[mfilename,'(''Quit'');']);
    set(data.handles.VariableDisplayed,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.Edges,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.PlotQuantity,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.PolarGrad,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.Thick,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.Black,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.AutoScale,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.Alinesize,'callback',[mfilename,'(''VariableDisplayed'');']);
    set(data.handles.Movie,'callback',[mfilename,'(''Movie'');']);
    set(data.handles.StartNumber,'callback',[mfilename,'(''StartNumber'',1);']);
    set(data.handles.EndNumber,'callback',[mfilename,'(''StartNumber'',2);']);
    set(data.handles.figMain,'CloseRequestFcn',[mfilename,'(''Quit'',0);']);
    set(data.handles.xyz,'String','');
    set(data.handles.Alinesize,'Value',1);
    data.startdirectory=pwd;
    data.current_figure=zeros(3,1);
    data.movieflags=zeros(3,1);
    data.playing=false;
    set(data.handles.Edges,'value',0);
    set(data.handles.EndNumber,'String',num2str(1));
    set(data.handles.StartNumber,'String',num2str(1));
    data.StartNumber=1;
    data.EndNumber=1;
    data=Project(data,true);
end

function data=viewmesh(data,option,projectnumber)
    if nargin<2
        option=0;
    end
    if nargin<3
        projectnumber=1;
    end
    if 0 %data.number==0
        data.m(projectnumber) = loadmesh_anyfile( [], fullfile(data.mdir{projectnumber},data.meshfile));
        data.m(projectnumber) = leaf_plot( data.m(projectnumber),...
            'morphogen',data.VariableDisplayed_i);
        data.handles.plotfigure=gcf;
        data.handles.plotaxis=gca;
        cameratoolbar
        string=[num2str(data.number),' (',num2str(data.m(projectnumber).globalProps.timestep),')'];
        set(data.handles.Time,'string',string);
        set(data.handles.Progress,'Value',data.number/max(data.meshnumbers));
    else
        loadedfile=[data.meshroot{projectnumber},makestagesuffixf(data.number),'.mat'];
        data.meshfile=fullfile(data.mdir{projectnumber},loadedfile);
        if exist(data.meshfile)==2 && option==0
            data.m(projectnumber).globalProps.mov = [];
            data.m(projectnumber).globalProps.makemovie = 0;
            data.m(projectnumber).globalProps.moviefile = '';
            data.m(projectnumber).globalProps.allowsave = 1;
            temp_plotdefaults    =data.m(projectnumber).plotdefaults;
            data.m(projectnumber).plothandles = struct();
            verbose=true;
            if projectnumber==1
%                 data.m(projectnumber)=loadmesh_anyfile( [], data.meshfile);
                temp_m=loadmesh_anyfile( [], data.meshfile);
                % ensure that fields are compatible
%                 temp_m=ReconcileFields(temp_m,data.m(projectnumber),verbose);
%                 data.m(projectnumber)=ReconcileFields(data.m(projectnumber),temp_m,verbose);
                data.m = temp_m;                
            else
                temp_m=loadmesh_anyfile( [], data.meshfile);
                % ensure that fields are compatible
                % new structure picks up fields from the old ones
%                 temp_m=ReconcileFields(temp_m,data.m(projectnumber-1),verbose);
                % and the old ones pick up fields from the new structure
%                 for temp_m_i=1:projectnumber-1
%                     data.m(temp_m_i)=ReconcileFields(data.m(temp_m_i),temp_m,verbose);
%                 end
                data.m(projectnumber) = temp_m;                
            end
            data.m(projectnumber).plotdefaults=temp_plotdefaults;
            data.m(projectnumber).plotdefaults.bioApointsize=0;
            if projectnumber==1 && isempty(get(data.handles.xyz,'String'))
                scale=max(abs(data.m(projectnumber).plotdefaults.axisRange));
                set(data.handles.xyz,'String',num2str(scale));
            end
% Commented out by RK as 'plotmode' no longer exists.  This
% appears to be a debugging statement: why do you want to dump
% the entire plot options here?
%             if ~isfield(data.m(projectnumber).plotdefaults,'plotmode')
%                 data.m(projectnumber).plotdefaults
%             end
        else
            set(data.handles.Play,'value',0);
        end
        if length(data.m)
            s=get(data.handles.PlotQuantity,'string');
            v=get(data.handles.PlotQuantity,'value')
            PlotQuantity=s{v};
            if data.current_figure(projectnumber)~=0
                figure(data.current_figure(projectnumber))
            end
            data.VariableDisplayed_i=get(data.handles.VariableDisplayed,'value');
            data.VariableDisplayed_n=get(data.handles.VariableDisplayed,'string');
            if data.current_figure(projectnumber)~=0 && ishandle(data.current_figure(projectnumber))
                figure(data.current_figure(projectnumber));
            end
% Commented out by RK as 'plotmode' no longer exists.  This
% appears to be a debugging statement: why do you want to dump
% the entire plot options here?
%             if ~isfield(data.m(projectnumber).plotdefaults,'plotmode')
%                 data.m(projectnumber).plotdefaults
%             end
            data.m(projectnumber).plotdefaults.bioAlinesize=get(data.handles.Alinesize,'Value');
            data.m(projectnumber).plotdefaults.bioApointsize=0;
            data.m(projectnumber).plotdefaults.thick=get(data.handles.Thick,'Value');
            data.m(projectnumber).plotdefaults.doClip=get(data.handles.doClip,'Value');
            if get(data.handles.Black,'Value')==0
                data.m(projectnumber).plotdefaults.bgcolor=[1,1,1];
            else
                data.m(projectnumber).plotdefaults.bgcolor=[0,0,0];
            end
            data.m(projectnumber).plotdefaults.autoScale=get(data.handles.AutoScale,'Value');
            scale=str2num(get(data.handles.xyz,'String'));
            data.m(projectnumber).plotdefaults.axisRange=[-scale, scale, -scale, scale, -scale, scale];
            q = leaf_plot( data.m(projectnumber),...
                'outputquantity', lower(PlotQuantity),...%'morphogen',...
                'morphogen',data.VariableDisplayed_i,...
                'drawedges', 2*get(data.handles.Edges,'value'),... % 0 or 2
                'drawgradients', get(data.handles.PolarGrad,'value')==1);%... % Boolean
            if length(data.m)>1
                data.m(projectnumber)=q;
            else
                data.m(projectnumber)=q;
            end
% Commented out by RK as 'plotmode' no longer exists.  This
% appears to be a debugging statement: why do you want to dump
% the entire plot options here?
%             if ~isfield(data.m(projectnumber).plotdefaults,'plotmode')
%                 data.m(projectnumber).plotdefaults
%             end
            data.current_figure(projectnumber)=gcf;
            data.current_axes(projectnumber)=gca;
%             string=sprintf('Step %5d DT=%3.3f Time=%7.3f sec Nodes %d',data.number,...
            string=sprintf('DT=%3.3f Time=%7.3f sec Nodes %d',...
                data.m(projectnumber).globalProps.timestep,...
                data.m(projectnumber).globalDynamicProps.currenttime,...
                size(data.m(projectnumber).nodes,1));
            set(data.handles.Time,'string',string,'FontWeight','bold');
            set(data.handles.Progress,'Value',data.number/max(data.meshnumbers));
            if get(data.handles.Black,'Value')~=0
                set(gca,'color','k');
            else
                set(gca,'Color','w');
            end
            drawnow
            if data.movieflags(projectnumber) &&  data.playing;
                if data.CurrentNumber<=data.EndNumber(projectnumber)
                    F=mygetframe(gca);
                    data.mov{projectnumber} = addframe(data.mov{projectnumber},F);
                else
                    if ishandle(data.mov{projectnumber})
                        data.mov{projectnumber}=close(data.mov{projectnumber});
                    end
                end
            end
            set(gcf,'name',loadedfile);
        end
    end
%     data=link_all_the_figures(data);
end

function data=Play(data,options)
    k=0;
%     v=get(data.handles.Progress,'Value');
%     k=ceil(v*max(data.meshnumbers));
    v=get(data.handles.Progress,'Value');
    n=ceil(v*max(data.meshnumbers));
    dist=(data.meshnumbers-n).^2;
    [dists,ii]=min(dist);
    k=ii;
    c=get(data.handles.Progress,'BackgroundColor');
    set(data.handles.Progress,'BackgroundColor','r');
    data.playing=true;
    while get(data.handles.Play,'value') ...
            && k<length(data.meshnumbers)
        k=k+1;
        data.number=data.meshnumbers(k);
        data.CurrentNumber=k;
        for i=1:length(data.m)
            data=viewmesh(data,0,i);
        end
    end
    if get(data.handles.Movie,'Value')
        set(data.handles.Movie,'Value',0);
        data=Movie(data,0);
    end
    data.playing=false;
    set(data.handles.Play,'value',0);
    set(data.handles.Progress,'BackgroundColor',c);
end

function data=Project(data,projectnumber)
    if nargin<2
        projectnumber=false;
    end
    if ~isfield(data,'mdir')
        [file,tempdir]=uigetfile(pwd,'Select a file in the project:movie:experiment_id directory');
        data.mdir{projectnumber}=tempdir;
    else
        pathname=fileparts(fileparts(data.mdir{1}));
        [file,data.mdir{projectnumber}]=uigetfile(pathname,'Select a file in the project:movie:experiment_id directory');
    end
    if file==0
        return;
    end
    commanddir1=data.mdir{projectnumber};
    commanddir2=fileparts(commanddir1(1:end-1));
    DD=dir(fullfile(commanddir2,'CommandLine.txt'));
%     DD=dir(fullfile(data.mdir{projectnumber},'CommandLine.txt'));
    if ~isempty(DD)
        h=fopen(fullfile(commanddir2,'CommandLine.txt'),'r');
        string=fgetl(h);
        fclose(h);
        if exist(fullfile(commanddir2,'ParameterSettings.txt'))==2
            h=fopen(fullfile(commanddir2,'ParameterSettings.txt'),'r');
            string=[string, ') [', fgetl(h),']'];
            fclose(h);
        end
        if projectnumber==1
            set(data.handles.Details1,'string',string,'FontWeight','bold');
        elseif projectnumber==2
            set(data.handles.Details2,'string',string,'FontWeight','bold','Visible','off');
        elseif projectnumber==3
            set(data.handles.Details3,'string',string,'FontWeight','bold','Visible','off');
        end
    end
    [data.pdir,proj] = fileparts(fileparts(data.mdir{projectnumber}));
    if projectnumber==1
        set(data.handles.Project1,'string',proj);
    elseif projectnumber==2
        set(data.handles.Project2,'string',proj);
    elseif projectnumber==3
        set(data.handles.Project3,'string',proj);
    end
  % data.meshroot(projectnumber)={file(1:strfind(file,gMISC_GLOBALS.stageprefix)-1)};%{file(1:end-8)};
  % Extract the part of the filename before the '_sNNNNdNNN.mat' suffix.
    data.meshroot(projectnumber) = ...
        {regexprep( file, ['^(.*)', gMISC_GLOBALS.stageprefix, '[0-9]+(d[0-9]*)?\.mat$'], '$1' )};
    fnstr=fullfile(data.mdir{projectnumber},['*', gMISC_GLOBALS.stageprefix, '*.mat']);
    DD=dir(fnstr);
    k=0;
    nn=[];
    for i=1:length(DD)
        nm=DD(i).name;
%         if isempty(findstr(nm,'0000000'))
            k=k+1;
            x = regexprep( nm, ['^.*', gMISC_GLOBALS.stageprefix, '([0-9]*(d[0-9]*)?)\.mat$'], '$1' );
            x = regexprep( x, 'd', '.' );
            nn(k) = str2num(x);
%             nm2=nm(strfind(nm,gMISC_GLOBALS.stageprefix)+2:end-4);
%             nm2(strfind(nm2,'d'))='.';
%             nn(k)=str2num(nm2);%nm(end-7:end-4));
%         end
    end
    data.meshnumbers=nn;
    data.number=nn(end);
    fnstr=fullfile(data.mdir{projectnumber},nm);
%     if ~isfield(data,'m')
        data.m(projectnumber) = loadmesh_anyfile( [], fnstr);
        %     data.number=0; % use zero to flag that system is being initialised
        ListOfMorphogens=fieldnames(data.m(projectnumber).mgenNameToIndex);
        data.number=nn(end);
        if ~iscell(get(data.handles.VariableDisplayed,'string'))
            set(data.handles.VariableDisplayed,'string',ListOfMorphogens,'value',1);
        end
        set(data.handles.Progress,'Value',1);
        set(data.handles.EndNumber,'String',num2str(max(data.meshnumbers)));
        data.EndNumber(projectnumber)=max(data.meshnumbers);
        if data.current_figure(projectnumber)~=0
            pos=get(data.current_figure(projectnumber),'Position');
            close(data.current_figure(projectnumber));
        end
        %     data=VariableDisplayed(data,0,projectnumber);
        %     data.current_figure(projectnumber)=gcf;
        %     data.current_axes(projectnumber)=gca;
%     end
    data=viewmesh(data,0,projectnumber);    
    %     data=VariableDisplayed(data,0,projectnumber);
        if exist('pos')
            set(data.current_figure(projectnumber),'Position',pos);
        end
    figno=gcf;
    if projectnumber==1
        set(data.handles.Project1,'string',['Fig.',num2str(figno),': ',proj]);
    elseif projectnumber==2
        set(data.handles.Project2,'string',['Fig.',num2str(figno),': ',proj]);
    elseif projectnumber==3
        set(data.handles.Project3,'string',['Fig.',num2str(figno),': ',proj]);
    end
end

function data=link_all_the_figures(data)
    list=[];
    for i=1:length(data.current_axes)
        if ishandle(data.current_axes(i))
            list(end+1)=data.current_axes(i);
        end
    end
   % set(gcf,'MenuBar','figure');
   list=list(:)';
    if length(data.current_axes)>1
        linkaxes(list);
        data.hlink=linkprop(list,{'CameraPosition','CameraUpVector'});
    else
        cameratoolbar(gcf);
    end
end

function data=Progress(data,options)
    v=get(data.handles.Progress,'Value');
    n=ceil(v*max(data.meshnumbers));
    dist=(data.meshnumbers-n).^2;
    [dists,ii]=min(dist);
    data.number=data.meshnumbers(ii);
    disp(sprintf('frame %f',data.number))
    set(data.handles.StartNumber,'String',num2str(data.number));
    if get(data.handles.Play,'value')==0
        for projectnumber=1:length(data.m)
            data=viewmesh(data,0,projectnumber);
        end
    end
end

function data=VariableDisplayed(data,flag,projectnumber)
    for projectnumber=1:length(data.m)
        data=viewmesh(data,flag,projectnumber);
    end
end
