function [argument_list,subfileArgList]=GFtboxCommand(varargin) %ProjectName,Directory,number_steps,stepsize)
    %function GFtboxCommand(ProjectName,Directory,number_steps,stepsize)
    %
    %For running growth models in batch mode. Results are filed in a
    %subdirectory
    %ProjectName:Movies:experimentname:meshes   stores full mesh .mat files
    %
    %At the end of each stage/step the mesh is stored in a separate file.
    %In addition, the command line used to invoke GFtboxCommand
    %and a copy of the interaction function are saved in the experimentname directory.
    %
    %The meshes can be viewed using VMS (View Mesh System)
    %or VMSReport. They can also be moved into the project directory where
    %GFtbox will recognise them as 'Stages'.
    %
    %When running jobs on the cluster, GFtboxCommand will automatically 
    %log the fact with a program called ClusterControl. 
    %
    %To run on the cluster you must have 
    %      1) Your (UEA) Cluster username 
    %      2) Your (UEA) Cluster username Password
    %      3) The (UEA) Cluster name, e.g. eslogin.uea.ac.uk
    %      4) Path to your filespace within the cluster (excluding your
    %      name)
    %             (e.g 'gpfs/tmpscratch01/sys/')
    %             (or  'gpfs/sys/')
    %and the following programs installed on your computer
    %      4) WinSCP
    %      5) Putty
    %      6) pscp
    %      7) plink
    %these can be installed on a Mac from the Terminal using MacPorts (e.g. sudo port install putty)
    %      
    %There is an option to output a set of shell scripts (.sh files) that
    %enable the results to be computed on a cluster in parallel.
    %This is invoked by the 'Use','Cluster', command pair.
    %
    %Usage, all argument pairs are optional
    %GFtboxCommand(...
    %    'State','RunNow',...   % status of command, options are:
    %                           % default is 'RunNow', i.e. when 'State' is not used then
    %                           %     system is run on the current computer
    %                           %
    %                           % The following only apply to running on the cluster
    %                           % 'Run'     create and run a single command
    %                           %
    %                           % 'Begin'   to begin assembling a batch of commands
    %                           % 'Add'     to add to the existing batch
    %                           % 'Start'   to add to and then run batch of commands
    %                           %
    %    'Use','Cluster'        % use the cluster if it is accessible
    %                           % default is not to use the cluster
    %    'Name','IMotif2a',...  % project to be run
    %    'Path', pwd,...         % present working directory
    %
    % either specify the modelling time in steps or stages
    %    'dt',[0.1, 0.01,...]   % default value is set by the project itself
    %                           % this command is used to change dt
    % Using time
    %    'Stages',[80,120,...]  % default, output at every iteration
    %                           % list of stages (in hours) that are to be
    %                           % output.
    %                           % Overrides Stepsize.
    %                           % Overrides total_time, i.e. total_time=max(stages)
    %                           % stage filenumbers are time - consistent
    %                           % GFtbox
    %    'Total_time',tvalue,...% i.e. number_steps =  tvalue/dt, also see Stages
    %    
    % Using steps
    %    'N_steps', 20,...      % number of steps to be run
    %    'Stepsize',1,...       % interval between steps that are recorded
    %                           % overridden by Stages
    %                           % defaut 1
    %
    %    'Flatten',true or false,...% The mesh will be dissected, (edges to be dissected must
    %                           % be entered into the mesh using the GUI).
    %                           % The mesh will then be dissected,exploded and flattened.
    %                           % after the last step (or stage).
    %                           % it can take a long time.
    %    'Mut_list',{'cych','div'},...% cell array of morphogen names to be set
    %                           %  to zero activity
    %                           % default is 'wildtype' i.e. all morph_a=1
    %    'param',[indices]      % it is possible to automatically
    %                           %create set of experiments that systematically vary parameters
    %                           %of the model. These have to be setup inside the interaction
    %                           %function. e.g. the parameter called 'param' is set up with
    %                           %m.userdata.ranges.param.range=[0, 1.2, 5, 10]
    %                           %m.userdata.ranges.param.index=3; % this sets the default, i.e.
    %                           %param=5. By setting an invalid index using
    %                           %on the GFtboxCommand command line, i.e. 10
    %                           in this case, each value of the parameter
    %                           will be used. If two range parameters are
    %                           set then the combinations will be run as
    %                           separate jobs, i.e. if there were 4 and 5
    %                           alternatives respectively, 20 experiments
    %                           would be created.
    %
    %                           %Sensitivity testing
    %                           %m.userdata.ranges.sensitivity.range
    %                           %m.userdata.ranges.sensitivity.index
    %                           this is a reserved userdatat ranges field
    %                           that is recognised by GFtboxCommand. In
    %                           To perform a sensitivity test, there must
    %                           be set of obj files in a directory called
    %                           objs. One obj per time step.
    %                           It is differences between the test sensitivity values
    %                           and these obj files that will be computed.
    %                           that will 
    %    'Subdiv',false,        % default, false,
    %                           % true then FE model will be subdivided
    %                           % to four times the number of vertices at the beginning
    %                           %
    %    'ExpID',name,          % not for general use. Used internally by GFtboxCommand
    %                           % to label experiments as they are run on
    %                           % the cluster
    % e.g.
    %
    %
    %Run on personal computer
    %GFtboxCommand('Name','HoleMotif4','Stages',[5,15,25]);
    %Run on cluster, exploring two range parameters
    %GFtboxCommand('Name','HoleMotif4','Stages',[5,15,25],'wound_production',3,'growth',1,'Use','Cluster');
    %
    %Run on cluster, exploring the effect of different dt and number of
    %vertices
    %GFtboxCommand('State','Begin','Name','NewFig4','dt',1.0,'Stages',[2,5,7],'Use','Cluster');
    %GFtboxCommand('State','Add','Name','NewFig4','dt',0.1,'Stages',[2,5,7],'Use','Cluster');
    %GFtboxCommand('State','Add','Name','NewFig4','dt',1.0,'Subdiv',true,Stages',[2,5,7],'Use','Cluster');
    %GFtboxCommand('State','Start','Name','NewFig4','dt',0.1,'Subdiv',true,'Stages',[2,5,7],'Use','Cluster');
    %
    %Try a range of dts, run for two separate models, for all combinations
    %of the specified mutations
    %GFtboxCommand('State','Start','Name','NewFig4','dt',[0.1 0.5 1] ,'Subdiv',true,'Stages',[2,5,7],'Use','Cluster','modelname',[1 2],'Mut_List',{'mutation1','mutation1'});
    %
    %J.Andrew Bangham, 2008
    %Jacob Newman, 2011
    %     if isunix
    %         % then assume it is running on the cluster and the paths are to be
    %         % setup relative to the current directory
    %         pathnames=fullfile(pwd,'DArT_Toolshed');
    %         addpath(genpath(pathnames));
    %     end
    %
    %

    MAKEPNG=false;
    args='';
    if rem(length(varargin),2)
        error('GFtboxCommand: arguments should be in pairs');
    end
    ProjectName='IMotif2a';
    if ispc
        result=system('plink.exe');
        if result~=1 %exist('plink.exe')~=0
            error('need to have plink.exe and pscp.exe on the system path, e.g. in C:\Program Files (x86)\Windows Live\Shared')
            %         if exist('C:\LocalStorage')==7
            %          path(path,'C:\LocalStorage');
            %         end
        end
    end
    Directory=pwd;
    number_steps=20;
    stepsize=1;
    dt=[];
    mutation_list={};
    subdivision_step_nums=false;
    total_time=[];
    range_names={};
    range_indexes={};
    flatten=[];
    runnow='RUNNOW';% 'RUN';
    stages=[];
    usecluster=false;
    ClearPasswords=false;
    makepng=false;
    makemesh=true;
    mov_VideoName='';
    disp(sprintf('Start - %s',datestr(now)));
    starttime=clock;
    short=true;
    sensitivityflag=false;
    for i=1:2:length(varargin)
        arg=upper(varargin{i});
        opt=varargin{i+1};
        switch arg
            % LOGISTICS
            case 'USE'
                usecluster=true;
            case 'PATH'
                Directory=opt;
            case 'NAME'
                ProjectName=opt;

            case 'PNG'
                makepng=opt;
            case 'MESH'
                makemesh=opt;
            case 'EXPID'
                % the following suffix is a daft legacy issue that arose when the the specification 
                % of stage filenames changed 
                if strcmp(opt(end-3:end),'_000')
                    mov_VideoName=opt;
                else
                    mov_VideoName=[opt,'_000']; % internal use (names have to be generated upfront
                    % to ensure no clashes due to concurrent processes
                end
            case 'CLEARPASSWORDS'
                ClearPasswords=true;

            case 'STATE'
                runnow=upper(opt);
                if isempty(strfind({'RUN','BEGIN','ADD','START','RUNNOW','RUNLATER'},runnow))
                    error('''State'', parameter should be one of ''RUN'',''BEGIN'',''ADD'',''START'',''RUNNOW'' ');
                end
                % NUMERICAL
                % either specify simulation duration using steps
            case 'N_STEPS'
                number_steps=opt;
            case 'STEPSZ'
                stepsize=opt;
                % or specify simulation duration using time
            case 'TOTAL_TIME',
                total_time=opt;
            case 'STAGES'
                stages=opt;
                % the relation between time and steps is
            case 'DT'
                dt=opt;
            case 'SHORT'
                if islogical(opt)
                    short=opt;
                elseif opt==1
                    short=true;
                else
                    short=false;
                end

            case 'SUBDIV'
                subdivision_step_nums=opt==true;
            case 'FLATTEN'
%                 flatten.start=opt(1);
                flatten=opt(1);
                
                % BIOLOGY
            case 'MUT_LIST'
                mutation_list=opt;
            otherwise % it is assumed that it is an interaction function range control
                %error(sprintf('GFtboxCommand: argument %s not recognised',arg))
                %it must be the name of an interaction function value that
                %is expected to be in a range,
                range_names{end+1}=varargin{i};
                if strcmpi(range_names{end},'sensitivity')
                    sensitivityflag=true;
                end
                %and the opt should EITHER
                %        be one valid index into the range specified in
                %        the interaction function
                %OR      be an index larger than the max possible and this
                %        indicates that all the values in the range are
                %        used
                range_indexes{end+1}=opt(:)';
        end
    end

    if exist(fullfile(Directory,ProjectName))~=7
        error(sprintf('SilentGFtool: cannot find %s in %s',...
            ProjectName,Directory));
    else % check that font case is correct (it matters on Unix)
        DD=dir(Directory);
        for i=1:length(DD)
            if strcmpi(DD(i).name,ProjectName)
                if ~strcmp(DD(i).name,ProjectName)
                    error(sprintf(...
                        '%s is not the same as the %s in the command, case matters on the unix cluster',...
                        DD(i).name,ProjectName));
                end
            end
        end
    end
    if usecluster 
        if strcmp(runnow,'RUNNOW')
            disp('This experiment will be submitted as a single job to the cluster');
            runnow='RUN';
        end
        %m = leaf_reload( m, 'restart' , 'rewrite', false);
        if strcmp(runnow,'RUN') || strcmp(runnow,'BEGIN') || exist('batchnumber.txt')~=2
            batchnumberstr=sprintf('%04d',floor(1000*rand));
            hbatch=fopen('batchnumber.txt','w');
            fprintf(hbatch,'%s',batchnumberstr);
            fclose(hbatch);    
        else
            hbatch=fopen('batchnumber.txt','r');
            batchnumberstr=fgetl(hbatch);
            fclose(hbatch);
        end
        % If the job is starting, then the batchnumber is no longer needed
        if strcmp(runnow,'START')
            delete('batchnumber.txt');
        end
        filename=sprintf('%s%s.sh',ProjectName,batchnumberstr);
        [errors,objflag]=Initial_SendToCluster(ProjectName,Directory,filename,ClearPasswords);
        if ~objflag && sensitivityflag
            error('A sensitivity test has been asked for, but there are no reference obj files')
        end
    end
    current_working_directory=pwd;
    cd(Directory);
    if exist(ProjectName)~=7
        error(sprintf('Unable to locate project %s in %s',ProjectName,Directory));
    end
    m = leaf_loadmodel( [], ProjectName, Directory, 'rewrite', false);
    m = leaf_setproperty( m, 'staticreadonly', true, 'allowInteraction',true ); % prevent cross talk when on cluster
    % do a sanity check on morphogen names if mutants have been enabled
    if ~isempty(mutation_list)
        for i=1:length(mutation_list)
            indx = FindMorphogenIndex( m, mutation_list{i}, 'Could not find morphogen');
            if isempty(indx)
                disp(sprintf('''MUT_LIST'',{''%s''}',mutation_list{i}))
                error('Mutant morphogen name is not valid for this model');
            end
        end
    end
    m = leaf_reload( m, 'restart' , 'rewrite', false);
    m = leaf_setproperty( m, 'staticreadonly', true, 'allowInteraction',true ); % prevent cross talk when on cluster
    % initialise the model
    % so that we can extract all range variables
    [m,OK] = leaf_dointeraction( m );
    if ~OK
        error('Problem with the interaction function');
    end
    % Decide on number of steps
    dt_list = [];
    if isempty(dt)
        dt=m.globalProps.timestep;
    elseif length(dt) > 1
        dt_list = dt;
        dt = dt(1);
    end
    timestep=dt; %Changed 25/05/2011

    % time overrides steps (steps are already set by default)
    if ~isempty(total_time)
        number_steps=round(total_time/dt);
    else
        total_time=number_steps*dt;
    end
    if ~isempty(stages)
        number_steps=round(max(stages)/dt);
        total_time=max(stages);
    else
        stages=0:stepsize*dt:total_time;
    end
    disp(sprintf('Computing %d intervals of %5.3f with outputs at',floor(total_time/dt),dt))
    disp(sprintf('%3.3f ',stages));
    disp('');
    %  %  m = leaf_iterate( m, 1, 'plot', 0 );
    %     set(gcf,'Visible','off');
    %     drawnow
    if isfield(m.userdata,'ranges')
        m_range_names=fieldnames(m.userdata.ranges);
        disp('The interaction function has ranges for variables')
        for i=1:length(m_range_names)
            disp(sprintf('%s ',m_range_names{i}));
        end
    else
        m_range_names='';
    end
    % Now have to loop through all the parameter ranges
    % checking each string in range_names against the list in ranges
    names={};
    for i=1:length(range_names)
        rangename_recognised=false;
        name=range_names{i};
        for j=1:length(m_range_names)
            if strcmpi(m_range_names{j},name)
                rangename_recognised=true;
                mfile_index=1:length(m.userdata.ranges.(m_range_names{j}).range);
                inters=intersect(range_indexes{i},mfile_index);
                if isempty(inters)%range_indexes{i}==0 % then loop through all
                    range{i}=1:length(m.userdata.ranges.(m_range_names{j}).range);
                    names{i}=name;
                elseif ~isempty(inters)
                    range{i}=inters;
                    names{i}=name;
                else
                    ss=sprintf('%f ',mfile_range);
                    disp(sprintf('ignore argument %s index [%d] interaction function range [%s]',...
                        name,range_indexes{i},ss));
                end
            end
        end
        if  ~rangename_recognised
            error(sprintf('%s range name not recognised',name));
        end
    end
    inputparameters=varargin;
    varin='';
    i=1;
    while i<=length(varargin)
        s=varargin{i};
        % must ignore when running the batch
        if strcmpi(s,'STATE') | strcmpi(s,'USE') | strcmpi(s,'DIR')
            i=i+1;
        else
            if ischar(s)
                varin=[varin,'''',s,''','];
            elseif iscell(s)
                varin=[varin,'{'];
                for ii=1:length(s)
                    varin=[varin,'''',s{ii},''','];
                end
                varin=[varin(1:end-1),'},'];
            else
                varin=[varin,'[',num2str(s),']',','];
            end
        end
        i=i+1;
    end
    varin(end)=' ';

    % Everything is now ready
    done_computing=false;
    if ~usecluster
        m = leaf_loadmodel( [], ProjectName, Directory, 'rewrite', false);
        m = leaf_setproperty( m, 'staticreadonly', true , 'allowInteraction',true); % prevent cross talk when on cluster
        m = leaf_reload( m, 'restart' , 'rewrite', false);
        m = leaf_setproperty( m, 'staticreadonly', true, 'allowInteraction',true ); % prevent cross talk when on cluster
        m = leaf_iterate( m, 1, 'plot', 0 ); % step past initialisation BECAUSE we want to change mutant number etc. from GFtboxCommand 
%         m = leaf_dointeraction( m ); % this would be safer were it not
%         for our having to handle all the mutantnumbers and modelnumbers
    end
    
    if ~isempty(mutation_list)
        % Remove the mut_list param from the command line args
        % These will be added later.
        indstart=(regexpi(varin,['''','mut_list','''','\s*',',{'],'start'));
        indend=indstart+regexpi(varin(indstart:end),['}'],'start');
        head=varin(1:indstart(1)-2);
        tail=varin(indend(1):end);
        varin=strtrim([head,tail]);
    end
    
    if ~isempty(dt_list)
        % Remove the dt param from the command line args
        % These will be added later.
        indstart=(regexpi(varin,['''','dt','''','\s*',',['],'start'));
        indend=indstart+regexpi(varin(indstart:end),[']'],'start');
        head=varin(1:indstart(1)-2);
        tail=varin(indend(1):end);
        varin=strtrim([head,tail]);
    end
    
    if ~isempty(names)
        % Remove the range params from the command line args
        % These will be added later.
        for i = 1:length(names)
            indstart=(regexpi(varin,['''',names{i},'''','\s*',',['],'start'));
            indend=indstart+regexpi(varin(indstart:end),[']'],'start');
            head=varin(1:indstart(1)-2);
            tail=varin(indend(1):end);
            varin=strtrim([head,tail]);
        end
    end
    
    % If there are no range variables, mutations or list of dts, then run
    % straight away.
    if isempty(names) && isempty(mutation_list) && isempty(dt_list) % no range variables
        if usecluster  %~strcmp(runnow,'RUN') %runnow>0
            disp(sprintf('Run later, no range variables specified, varin=%s %s',varin,mov_VideoName))
            temp_var=varin;
            argument_list{1}=temp_var;
        else
            disp(sprintf('Run here now, no range variables, varin=%s %s',varin,mov_VideoName))
            % no parameters to be replaced.
            [m,NewDir]=change_and_compute(m,timestep,Directory,ProjectName,...
                mutation_list,stepsize,number_steps-1,stages,...
                subdivision_step_nums,inputparameters,makepng,makemesh,mov_VideoName,flatten);
            done_computing=true;
        end
        fprintf(1,'\n');
    else
        % Changed 23/05/2011
        % If there are range variables specified.
        if ~isempty(names)
            % Generate all range permutations.
            disp(sprintf('Run later, range variables, varin=%s %s',varin,mov_VideoName))
            variants = range{1}';
            for i = 2:length(range);
                variants = repmat(variants,length(range{i}),1);
                variants(:,size(variants,2)+1) = reshape(repmat(range{i},size(variants,1)/length(range{i}),1),size(variants,1),1);
            end
            argument_list={};
            temp_var = varin;
        else
            % No ranges.
            disp(sprintf('Run later, mutations list, varin=%s %s',varin,mov_VideoName))
            variants = {};
            argument_list={varin};
        end
        if usecluster
            % Make parameters list for all permutations of ranges.
            for i = 1:size(variants,1)
                configTrace = '';
                newArg = temp_var;
                for n = 1:size(variants,2)
                    newArg=[newArg,',''',names{n},''',[',num2str(variants(i,n)),']'];
                    configTrace = [configTrace, ' {', names{n}, ' ', num2str(variants(i,n)),'}'];
                end
                fprintf(1,'>>>>>>>>>%s\n',configTrace);
                argument_list{end+1} = newArg;
            end
            % Prepare commands for all combinations of mutation_list, for
            % each range permuation.
            oldarglist = argument_list;
            for i = 1:length(mutation_list)
                combinations = nchoosek(1:length(mutation_list),i);
                for n = 1:size(combinations,1)
                    temp_mut_list = cellstr(strcat(' ''',char(mutation_list{combinations(n,:)}),''''));
                    temp_mut_list = strtrim(horzcat(temp_mut_list{:}));
                    mut_comb_string = sprintf(',''Mut_list'',{%s}',temp_mut_list');
                    argument_list = horzcat(argument_list,strcat(oldarglist',repmat(mut_comb_string,length(oldarglist),1))');
                    fprintf(1,'>>>>>>>>>%s\n',mut_comb_string(2:end));
                end
            end
            % Prepare commands for the list of dts.
            oldarglist = argument_list;
            if ~isempty(dt_list)
                argument_list = []; % We don't a set of commands which the default DT, if DT values have been specified
            end
            for i = 1:length(dt_list)
                    dt_string = sprintf(',''dt'',[%f]',dt_list(i));
                    argument_list = horzcat(argument_list,strcat(oldarglist',repmat(dt_string,length(oldarglist),1))');
                    fprintf(1,'>>>>>>>>>%s\n',dt_string(2:end));
            end
            fprintf(1,'\n');
        else
            % Run the range variable experiments locally.
            disp(sprintf('Run here now, range variables and/or mutations list, varin=%s %s',varin,mov_VideoName))
            % If there's no dt list, we want to use the default value or the single
            % specified value.
            if isempty(dt_list)
               dt_list = dt; 
            end
            if isempty(variants)
                variants = -1; % Value for variants if no range variables are specified 
            end
            % Go through the variants for the ranges first.
            for i = 1:size(variants,1)
                temp_var = varin;
                if variants(1) ~= -1
                    configTrace = '';
                    for n = 1:size(variants,2)
                        m.userdata.ranges.(names{n}).index=variants(i,n);
                        configTrace = [configTrace, ' {', names{n}, ' ', num2str(variants(i,n)),'}'];
                        temp_var = [temp_var, ',''', char(names{n}),''',[',num2str(variants(i,n)),']'];
                    end
                    fprintf(1,'>>>>>>>>>%s\n',configTrace);
                end
                % Go through the list of dts
                for nn = 1:length(dt_list)
                    % time overrides steps (steps are already set by
                    % default).
                    dt = dt_list(nn);
                    temp_var_dt = [temp_var, ',dt,[',num2str(dt),']'];
                    if ~isempty(total_time)
                        number_steps=round(total_time/dt);
                    else
                        total_time=number_steps*dt;
                    end
                    if ~isempty(stages)
                        number_steps=round(max(stages)/dt);
                        total_time=max(stages);
                    else
                        stages=0:stepsize*dt:total_time;
                    end
                    fprintf(1,'>>>>>>>>>dt,[%d]\n',dt_list(nn));
                    % Go through the mutation list, if there is one.
                    if isempty(mutation_list)
                        % Calc stepsize, number_steps etc
                        [m,NewDir]=change_and_compute(m,timestep,Directory,ProjectName,...
                            mutation_list,stepsize,number_steps,stages,...
                            subdivision_step_nums,temp_var_dt,makepng,makemesh,mov_VideoName,flatten);
                        done_computing=true;
                    else
                        % If there is a mutation list, calculate the
                        % combinations of the mutations. Start with no
                        % mutations.
                        temp_mut_list = {};
                        fprintf(1,'>>>>>>>>>No mutations\n');
                        [m,NewDir]=change_and_compute(m,timestep,Directory,ProjectName,...
                            temp_mut_list,stepsize,number_steps,stages,...
                            subdivision_step_nums,temp_var_dt,makepng,makemesh,mov_VideoName,flatten);
                        done_computing=true;
                        for ii = 1:length(mutation_list)
                            combinations = nchoosek(1:length(mutation_list),ii);
                            for n = 1:size(combinations,1)
                                temp_mut_list = cellstr(strcat(char(mutation_list{combinations(n,:)})))';
                                configTrace = cellstr(strcat(char(mutation_list{combinations(n,:)}),','))';
                                configTrace = horzcat(configTrace{:});
                                temp_var_mut = [temp_var_dt, ',Mut_list,{',configTrace(1:end-1),'}'];
                                fprintf(1,'>>>>>>>>>%s\n',configTrace(1:end-1));
                                [m,NewDir]=change_and_compute(m,timestep,Directory,ProjectName,...
                                    temp_mut_list,stepsize,number_steps,stages,...
                                    subdivision_step_nums,temp_var_mut,makepng,makemesh,mov_VideoName,flatten);
                                done_computing=true;
                            end
                        end
                    end
                end
            end
        end
    end
    fprintf(1,'\n');

    % run later
    if usecluster %isunix
%         disp('Run later on Cluster')
%         %             ping_cluster('ping escluster.uea.ac.uk');
%         if strcmp(runnow,'RUN') || strcmp(runnow,'BEGIN')
%             batchnumberstr=sprintf('%04d',floor(1000*rand));
%             hbatch=fopen('batchnumber.txt','w');
%             fprintf(hbatch,'%s',batchnumberstr);
%             fclose(hbatch);
%         else
%             hbatch=fopen('batchnumber.txt','r');
%             batchnumberstr=fgetl(hbatch);
%             fclose(hbatch);
%         end
        filename=sprintf('%s%s.sh',ProjectName,batchnumberstr);
        reportfilename=sprintf('%s%s.txt',ProjectName,batchnumberstr);
        NoExistingExp=0;
        %             if ~(strcmp(runnow,'RUN') || strcmp(runnow,'START'))
        if (strcmp(runnow,'ADD') || strcmp(runnow,'START'))
            % how many already in the list?
            if exist(filename)==2
                h=fopen(filename,'r');
                while ~feof(h)
                    s=fgets(h);
                    if strfind(s,'qsub')
                        NoExistingExp=NoExistingExp+1; % number of existing experiments
                    end
                end
                fclose(h);
            end
        end
        errors=false;
        for i=1:length(argument_list)
            ID=[ProjectName,batchnumberstr,'R',sprintf('%03d',(i+NoExistingExp)),'R','_000'];
%             ID(strfind(ID,' '))='-';
%             ID(strfind(ID,':'))='-';
            temp_argument_list=[argument_list{i},',''ExpID'',''',ID,''''];
            disp(temp_argument_list)
            %                 disp(argument_list{i})
            % send individual shell scripts to cluster, one per experiment (job)
            [errors_temp,subfilename]=...
                unixtemplate(NoExistingExp+i,temp_argument_list,usecluster,Directory,ProjectName,batchnumberstr,short);
            %                     unixtemplate(k+i,argument_list{i},usecluster,Directory,ProjectName);
            subfilenames{i}=subfilename;
            subfileArgList{i}=temp_argument_list;
            errors=errors | errors_temp;
        end
        % construct the shell script that will submit the jobs
        if strcmp(runnow,'RUN') || strcmp(runnow,'BEGIN')
            h_master=fopen(filename,'w');
            h_report=fopen(reportfilename,'w');
            fprintf(h_master,'#!/bin/csh\r\n');
        else
            h_master=fopen(filename,'a');
            h_report=fopen(reportfilename,'a');
        end
        errors=false;
        for i=1:length(argument_list)
            fprintf(h_master,'qsub %s\r\n',subfilenames{i});
            fprintf(h_report,'%s,\t,GFtboxCommand(%s);\n',subfilenames{i},subfileArgList{i});
        end
        fclose(h_master);
        fclose(h_report);
        % now run it
        if strcmp(runnow,'RUN') || strcmp(runnow,'START')
            if usecluster
                % now force a restart to ensure that the static file is
                % initialised
                errors=SendToCluster(ProjectName,Directory,filename,ClearPasswords);
            else
                RunSilent
            end
        end
    elseif  done_computing==false
        % run now from a batch file
        disp('Run later here')
        disp(sprintf('ProjectName=%s runnow=%s',ProjectName,runnow))
        filename=sprintf('%sRun.m',ProjectName);
        k=0;
        if strcmp(runnow,'RUN') || strcmp(runnow,'BEGIN')
            h_master=fopen(filename,'w');
        else % add and start
            % how many already in the list?
            if exist(filename)==2
                h=fopen(filename,'r');
                while ~feof(h)
                    s=fgets(h);
                    k=k+1;
                end
                fclose(h);
            end
            h_master=fopen(filename,'a');
        end
        for i=1:length(argument_list)
            disp(argument_list{i})
            mtemplate(i+k,argument_list{i},ProjectName);
            fprintf(h_master,'%sRun%-d\r\n',ProjectName,k+i);
        end
        fclose(h_master);
        if strcmp(runnow,'RUN') || strcmp(runnow,'START') %runnow==1 || runnow==4
            disp(sprintf(...
                'run the newly created Matlab script: %s first check it is present',filename))
            which(filename)
            feval(filename(1:end-2));
        end
    end
    cd(current_working_directory);
    disp(sprintf('Finish - %s',datestr(now)));
    disp(sprintf('Elapsed time = %d seconds',round(etime(clock,starttime))));
end

function  errors=SendToCluster(ProjectName,Directory,filename,ClearPasswords)
%     if ClearPasswords
%         clearpasswords
%     end
%     %ClusterSVN; %
% 
%     errors=false;
%     rs=which('RunSilent');
%     string0=sprintf('pscp -pw %s "%s" "%s@escluster.uea.ac.uk:/gpfs/%s%s/RunSilent.m"',...
%         getpassword,rs,getyourclustername,getyourclusterdirectory,getyourclustername);
%     errors=unixcluster(string0,'Copied RunSilent.m to current unix directory ','copy','RunSilent.m');
% 
%     % copy project to cluster (without movie directory)
%     string2=sprintf('plink -pw %s %s@escluster.uea.ac.uk cd /gpfs/%s%s/; mkdir %s;',...
%         getpassword,getyourclustername,getyourclusterdirectory,getyourclustername,ProjectName);
%     errors=unixcluster(string2,'Making project directory if needed ');
%     interactionpath=Directory;
%     interactionfunction=fullfile(ProjectName,[lower(ProjectName),'.m']);
%     unix_interactionfunction=[ProjectName,'/',lower(ProjectName),'.m'];
%     string0=sprintf('pscp -pw %s "%s%s%s" "%s@escluster.uea.ac.uk:/gpfs/%s%s/%s/%s.m"',...
%         getpassword,interactionpath,filesep,interactionfunction,getyourclustername,...
%         getyourclusterdirectory,getyourclustername,ProjectName,lower(ProjectName));
%     errors=unixcluster(string0,'Copied interaction function to ','copy',interactionfunction);
% 
%     string0=sprintf('pscp -pw %s "%s" "%s@escluster.uea.ac.uk:/gpfs/%s%s/%s/%s.mat"',...
%         getpassword,fullfile(Directory,ProjectName,[ProjectName,'.mat']),getyourclustername,...
%         getyourclusterdirectory,getyourclustername,ProjectName,ProjectName);
%     disp('The next copy can take a minute or so')
%     errors=unixcluster(string0,'Copied structure file to ','copy',[ProjectName,'.mat']);
% 
%     string0=sprintf('pscp -pw %s "%s" "%s@escluster.uea.ac.uk:/gpfs/%s%s/%s/%s_static.mat"',...
%         getpassword,fullfile(Directory,ProjectName,[ProjectName,'.mat']),getyourclustername,...
%         getyourclusterdirectory,getyourclustername,ProjectName,ProjectName);
%     disp('The next copy can take a minute or so')
%     errors=unixcluster(string0,'Copied static structure file to ','copy',[ProjectName,'.mat']);
% 
%     % create movies directory
%     string2=sprintf('plink -pw %s %s@escluster.uea.ac.uk cd /gpfs/%s%s/%s; mkdir movies',...
%         getpassword,getyourclustername,getyourclusterdirectory,getyourclustername,ProjectName);
%     errors=unixcluster(string2,'Making project movies directory if needed ');

    % copy master file to escluster
   
    string=sprintf('pscp -pw %s "%s" "%s@%s:/%s%s/%s"',...
        getpassword,fullfile(Directory,filename),getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
    errors=unixcluster(string,sprintf('Copied %ss to ',filename),'copy',filename);
    if ispc == 1
        string2=sprintf('plink -pw %s %s@%s cd /%s%s/; chmod +wx %s; dos2unix %s;',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename,filename);
    else
        string2=sprintf('plink -pw %s %s@%s chmod +wx /%s%s/%s',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
        errors=unixcluster(string2,'Make executable');
        string2=sprintf('plink -pw %s %s@%s dos2unix /%s%s/%s',getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
        errors=unixcluster(string2,'Make dos2unix');
    end
    errors=unixcluster(string2,'Make executable and dos2unix ');

    % and run the master shell script file which will run the others
    if ispc == 1
        string2=sprintf('plink -pw %s %s@%s cd /%s%s/; %s;',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
    else
        string2=sprintf('plink -pw %s %s@%s /%s%s/%s',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
    end
    errors=unixcluster(string2,sprintf('Start running %s ',filename));

    % AND create the retrieval file
    %create_retrieval_function(Directory,ProjectName,filename);

    % FINALLY report to ClusterControl
    % New 12/05/11 - ClusterBuffer is edited remotely
    ClusterControlAdd('Project',filename,'Directory',Directory);
    % Send the ClusterBuffer.txt file to cluster in case data is to be retrieved elsewhere
%     if ispc == 1
%         string=sprintf('pscp -pw %s "%s" "%s@%s:/gpfs/%s%s/ClusterBuffer.txt"',...
%         getpassword,fullfile(Directory,'ClusterBuffer.txt'),getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername);
%     else
%         string=sprintf('pscp -pw %s "%s" "%s@%s:/gpfs/%s%s/ClusterBuffer.txt"',...
%         password,fullfile(Directory,'ClusterBuffer.txt'),username,clustername,clusterDir,username); 
%     end
%     errors=unixcluster(string,sprintf('Copied %ss to ClusterBuffer.txt'),'copy','ClusterBuffer.txt');
    
    % AND remove the master shell script file
    newProjectName=[filename,'.done'];
    if ispc == 1
        [a,b]=system(sprintf('del %s',filename));
    else
        [a,b]=system(sprintf('rm %s',filename));
    end
    if a~=0
        disp(b);
    else
        disp(sprintf('deleted %s ',filename));
    end
%     % AND rename the master shell script file
%     newProjectName=[filename,'.done'];
%     [a,b]=system(sprintf('rename %s %s',filename,newProjectName));
%     if a~=0
%         disp(b);
%     else
%         disp(sprintf('renamed %s as %s',filename,newProjectName));
%     end
end


function  [errors,objflag]=Initial_SendToCluster(ProjectName,Directory,filename,ClearPasswords)
    if ClearPasswords
        clearpasswords
    end
    % If on Mac, install Putty for pscp and plink (via Mac Ports)
    % Also add
    % setenv('PATH','/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin:/usr/X11/bin:/opt/local/bin');
    %ClusterSVN; %
    objflag=false;
    errors=false;
    rs=which('RunSilent');
    string0=sprintf('pscp -pw %s "%s" "%s@%s:/%s%s/RunSilent.m"',getpassword,rs,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername);
    errors=unixcluster(string0,'Copied RunSilent.m to current unix directory ','copy','RunSilent.m');
    
    % copy project to cluster (without movie directory)
    if ispc == 1
    string2=sprintf('plink -pw %s %s@%s cd /%s%s/; mkdir %s;',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName);
    else
        string2=sprintf('plink -pw %s %s@%s mkdir /%s%s/%s',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName);
    end
    errors=unixcluster(string2,'Making project directory if needed ');
    interactionpath=Directory;
    interactionfunction=fullfile(ProjectName,[lower(ProjectName),'.m']);
    unix_interactionfunction=[ProjectName,'/',lower(ProjectName),'.m'];
    string0=sprintf('pscp -pw %s "%s%s%s" "%s@%s:/%s%s/%s/%s.m"',...
        getpassword,interactionpath,filesep,interactionfunction,getyourclustername,...
        getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName,lower(ProjectName));
    errors=unixcluster(string0,'Copied interaction function to ','copy',interactionfunction);
    if errors==true
        error('Failed to copy interaction function (could be that the name has - instead of _)');
    end
    string0=sprintf('pscp -pw %s "%s" "%s@%s:/%s%s/%s/%s.mat"',...
        getpassword,fullfile(Directory,ProjectName,[ProjectName,'.mat']),getyourclustername,...
        getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName,ProjectName);
    disp('The next copy can take a minute or so')
    errors=unixcluster(string0,'Copied structure file to ','copy',[ProjectName,'.mat']);
    string0=sprintf('pscp -pw %s "%s" "%s@%s:/%s%s/%s/%s_static.mat"',...
        getpassword,fullfile(Directory,ProjectName,[ProjectName,'_static.mat']),getyourclustername,...
        getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName,ProjectName);
    disp('The next copy can take a minute or so')
    errors=unixcluster(string0,'Copied static structure file to ','copy',[ProjectName,'.mat']);

    if exist(fullfile(Directory,ProjectName,'objs'))==7
        % copy project to cluster (without movie directory)
        if ispc == 1
            string2=sprintf('plink -pw %s %s@%s cd /%s%s/; ls %s/%s;',...
            getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName,'objs');
        else
            string2=sprintf('plink -pw %s %s@%s ls /%s%s/%s/%s',...
            getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName,'objs');
        end
        errors=unixcluster(string2,'Checking to see if obj directory already exists ');
        if errors
            % then make the directory and copy the obj files into place
            if ispc == 1
                string2=sprintf('plink -pw %s %s@%s cd /%s%s/; mkdir %s/%s;',...
                getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName,'objs');
            else
                string2=sprintf('plink -pw %s %s@%s mkdir /%s%s/%s/%s;',...
                getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName,'objs');
            end
            errors=unixcluster(string2,'Making project directory if needed ');
            D=dir(fullfile(Directory,ProjectName,'objs','*.obj'));
            for i=1:length(D)
                objflag=true;
                objfilename=D(i).name;
                string0=sprintf('pscp -pw %s "%s" "%s@%s:/%s%s/%s/objs/%s"',...
                    getpassword,fullfile(Directory,ProjectName,'objs',objfilename),getyourclustername,...
                    clustername,getyourclusterdirectory,getyourclustername,ProjectName,objfilename);
                disp(sprintf('Copying obj file %s',objfilename))
                errors=unixcluster(string0,'Copied obj file to ','copy',[ProjectName,'.mat']);
            end
        end
    end
    % create runs directory
    if ispc == 1
        string2=sprintf('plink -pw %s %s@%s cd /%s%s/%s; mkdir runs',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName);
    else
        string2=sprintf('plink -pw %s %s@%s mkdir /%s%s/%s/runs',...
        getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,ProjectName);

    end
    errors=unixcluster(string2,'Making project runs directory if needed ');
end

function [m,NewDir]=change_and_compute(m,timestep,Directory,ProjectName,...
        mutation_list,stepsize,number_steps,stages,...
        subdivision_step_nums,inputparameters,makepng,makemesh,mov_VideoName,flatten)
    %replace altered fields
    [mov_VideoName,NewDir] = setup_movie_folder_and_names(Directory,ProjectName,mov_VideoName);
    % HOUSEKEEPING, setup movie folder and
    keep_record_of_experiment(ProjectName,Directory,NewDir,inputparameters);
    m.globalProps.mov.VideoName=mov_VideoName;
    m.globalProps.timestep=timestep;
    m = leaf_setproperty( m, 'staticreadonly', true , 'allowInteraction',true ); % prevent cross talk when on cluster
    % Decide on mutations
    m = leaf_allowmutant( m, false ); % Warning, argument has to be logical, I copied a 1 from somewhere
    if ~isempty(mutation_list)
        m = leaf_allowmutant( m, true ); % Warning, argument has to be logical, I copied a 1 from somewhere
        m.mutantLevel=ones(size(m.mutantLevel)); % convert all to wildtype
        for i=1:length(mutation_list)
            indx = FindMorphogenIndex( m, mutation_list{i}, 'Could not find morphogen');
            m = leaf_mgen_modulate( m, 'morphogen', mutation_list{i}, 'mutant', 0 );
            disp(sprintf('Set %s to 0',mutation_list{i}))
        end
    else
        disp('Wild type')
    end
    m=compute_results(m,stepsize,number_steps,subdivision_step_nums,...
        stages,makepng,makemesh,NewDir,flatten);
end

function m=compute_results(m,stepsize,number_steps,subdivision_step_nums,stages,...
        makepng,makemesh,NewDir,flatten)
    
    % seed the random number using time and filename (to ensure different starting point for each
    % instance running on the cluster
    m = leaf_setproperty( m, 'staticreadonly', true, 'allowInteraction',true ); % prevent cross talk when on cluster
    seednumber=sum(10000*clock)+sum(double(NewDir));
    rand('twister',seednumber);
    randn('state',seednumber);
    m=leaf_record_mesh_frame(m,'RECORD','ON');
    if makemesh
        disp('save mesh')
        m=savemesh(m,NewDir);
    end
    if makepng
        disp('save png')
        m=savepng(m,NewDir)
    end

    if isempty(stages)
        last_step=0;
        % do two iterations before subdivision to allow any subdivisions in interaction function to
        % have completed
        m = leaf_iterate( m, 2, 'plot', 0 );
        if subdivision_step_nums %~isempty(subdivision_step_nums)
            m = leaf_dointeraction( m, 1 ); % to ensure that seams are already setup
            m = leaf_refineFEM( m, 1 );
        end
        % make quite sure that growth, diffusion and elasticity are switched on
        m.globalProps.growthEnabled=1;
        m.globalProps.diffusionEnabled=true;
        m.globalProps.plasticGrowth=1;
        for step=3:stepsize:number_steps
            realtime = m.globalDynamicProps.currenttime;
            last_step=step;
            % iterate stepsize steps without recording
            disp(sprintf('\n%s  number %d/%d dt=%f realtime=%f secs meshsize %d elements',...
                m.globalProps.mov.VideoName(1:end-4),step,round(number_steps/stepsize),m.globalProps.timestep,...
                m.globalDynamicProps.currenttime,size(m.nodes,1)))
            m = leaf_iterate( m, stepsize, 'plot', 0 );
            disp('save output');
            if makemesh
                disp('save mesh')
                m=savemesh(m,NewDir);
            end
            if makepng
                disp('save png')
                m=savepng(m,NewDir)
            end
        end
    else
        last_step=0;
        stepsize=1;
        realtime = m.globalDynamicProps.currenttime;
        step=1;
        % do two iterations before subdivision to allow any subdivisions in interaction function to
        % have completed
        m = leaf_iterate( m, 2, 'plot', 0 );
        if subdivision_step_nums %~isempty(subdivision_step_nums)
            m = leaf_dointeraction( m, 1 );% to ensure that seams are already setup
            m = leaf_refineFEM( m, 1 );
            m.globalProps.allowBentLongFEM=0;
            m.globalProps.allowSplitBentFEM=0;
            m.globalProps.allowBentLongFEM=0;
        end
        for i_stage=1:length(stages)
            while realtime<stages(i_stage)
                last_step=step;
                % iterate stepsize steps without recording
                disp(sprintf('\n%s  stage %d/%d realtime=%f dt=%f realtime=%f secs meshsize %d elements',...
                    m.globalProps.mov.VideoName(1:end-4),i_stage,length(stages),realtime,m.globalProps.timestep,...
                    m.globalDynamicProps.currenttime,size(m.nodes,1)))
                m = leaf_iterate( m, 1, 'plot', 0 );
                realtime = m.globalDynamicProps.currenttime;
            end
            disp('save output');
            if makemesh
                disp('save mesh')
                m=savemesh(m,NewDir);
            end
            if makepng
                disp('save png')
                m=savepng(m,NewDir)
            end
        end
    end
    if ~isempty(flatten)
        dt=m.globalProps.timestep;
        m.globalProps.timestep=m.globalProps.timestep/5;
        realtime = m.globalDynamicProps.currenttime;
        m_temp=m;
        m = leaf_iterate( m, 1, 'plot', 0 );
        realtime_after_first_small_step = m.globalDynamicProps.currenttime;
%         m.globalDynamicProps.currenttime=realtime + dt; % just pretending we have moved on
        m=leaf_dissect(m);
        m=leaf_flatten(m);
        disp(sprintf('---------------- first dissection NewDir=%s realtime=%f',NewDir,realtime_after_first_small_step));
        if makemesh
            disp('save mesh')
            m=savemesh(m,NewDir);
        end
        if makepng
            disp('save png')
            m=savepng(m,NewDir)
        end
%         m=m_temp;
%         m = leaf_iterate( m, 1, 'plot', 0 );
% %         m.globalDynamicProps.currenttime=realtime + 2*dt; % just pretending we have moved on
%         m=leaf_dissect(m);
%         m=leaf_flatten(m,1);
%         disp(sprintf('---------------- second dissection NewDir=%s realtime=%f',NewDir,realtime_after_first_small_step));
%         if makemesh
%             disp('save mesh')
%             m=savemesh(m,NewDir);
%         end
%         if makepng
%             disp('save png')
%             m=savepng(m,NewDir)
%         end
%         m=m_temp;
%         m = leaf_iterate( m, 1, 'plot', 0 );
% %         m.globalDynamicProps.currenttime=realtime + 2*dt; % just pretending we have moved on
%         m=leaf_dissect(m);
%         m=leaf_flatten(m,1);
%         disp(sprintf('---------------- third dissection NewDir=%s realtime=%f',NewDir,realtime_after_first_small_step));
%         if makemesh
%             disp('save mesh')
%             m=savemesh(m,NewDir);
%         end
%         if makepng
%             disp('save png')
%             m=savepng(m,NewDir)
%         end
%         number_stepsf=flatten/m.globalProps.timestep;
%         number_stages_viewed=3;
%         number_stepsf=round(number_stepsf/number_stages_viewed);
%         for i=1:number_stages_viewed
%             m = leaf_iterate( m, number_stepsf, 'plot', 0 );
%             disp('save dissected, exploded and flatted mesh');
%             if makemesh
%                 disp('save mesh')
%                 m=savemesh(m,NewDir);
%             end
%             if makepng
%                 disp('save png')
%                 m=savepng(m,NewDir)
%             end
%         end
    end
end


function [filename]=mtemplate(n,s,ProjectName)
    filename=sprintf('%sRun%-d.m',ProjectName,n);
    h=fopen(filename,'w');
    fprintf(h,'RunSilent(%s)\r\n',s);
    fclose(h);
end

function [errors,filename]=pctemplate(n,s,ProjectName)
    filename=sprintf('%sRun%-d.bat',ProjectName,n);
    h=fopen(filename,'w');
    fprintf(h,'matlab -nosplash -nodesktop -r "RunSilent(%s)"\r\n',s);
    fclose(h);
end

function [errors,filename]=unixtemplate(n,s,usecluster,Directory,ProjectName,batchnumberstr,short)
    if nargin<7
        short=true;
    end
    
    filename=sprintf('%s%sR%-d.sh',ProjectName,batchnumberstr,n);
    h=fopen(filename,'w');
    fprintf(h,'#!/bin/csh\r\n');
    if short==true
        fprintf(h,'#$ -q short.q\r\n');
    else
        fprintf(h,'#$ -q long.q\r\n');
    end
%     fprintf(h,'#$ -m e\r\n'); % mail me
%     fprintf(h,'#$ -q long.q\r\n');
    fprintf(h,'#$ -cwd\r\n');
    fprintf(h,'pwd\r\n');
    fprintf(h,'date\r\n');
    fprintf(h,'#$ -cwd\r\n');
    fprintf(h,'echo -- start time --\r\n');
    fprintf(h,'module add shared matlab/2008a\r\n');
%     fprintf(h,'#$ -M %s\r\n',getyouremail); % email on job completion
    fprintf(h,'matlab -nosplash -nodesktop -r "RunSilent(%s)"\r\n',s);
    fprintf(h,'echo -- end time --\r\n');
    fprintf(h,'date\r\n');
    fclose(h);
    if usecluster % need to send the newly created file to the cluster and set properties
        string=sprintf('pscp -pw %s "%s" "%s@%s:/%s%s/%s"',...
            getpassword,fullfile(Directory,filename),getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
        errors=unixcluster(string,sprintf('Copied experiment %s',filename),'move',filename);
        if ispc == 1
            string2=sprintf('plink -pw %s %s@%s cd /%s%s/; chmod +wx %s; dos2unix %s;',...
            getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename,filename);
            errors=unixcluster(string2,'Make executable and dos2unix ');
        else
            string2=sprintf('plink -pw %s %s@%s chmod +wx /%s%s/%s',...
            getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
            errors=unixcluster(string2,'Make executable');
            string2=sprintf('plink -pw %s %s@%s dos2unix /%s%s/%s',getpassword,getyourclustername,getTheClusterNameAddress,getyourclusterdirectory,getyourclustername,filename);
            errors=unixcluster(string2,'Make dos2unix');
        end
        disp('')
    end
end

function m=savemesh(m,NewDir);
    saveasbase = (m.globalProps.allowsave && (m.globalDynamicProps.currentIter==0));
    if saveasbase==0
        correct_VideoName=m.globalProps.mov.VideoName;

        default_plot_parameters.pictures=m.pictures;
        default_plot_parameters.uicontrols=m.plotdefaults.uicontrols;
        default_plot_parameters.invisibleplot=m.plotdefaults.invisibleplot;

        m=temporary_replacement_of_name(m,'meshes',NewDir);
        m.pictures=[];
        m.plotdefaults.invisibleplot=false;
        m.plotdefaults.uicontrols=true;
        disp(sprintf('About to save mesh'))
        m=leaf_record_mesh_frame(m);
        makestagesuffixf( m.globalDynamicProps.currenttime );
        disp(sprintf('written %s',[m.globalProps.modelname,makestagesuffixf( m.globalDynamicProps.currenttime )]))
        m.globalProps.mov.VideoName=correct_VideoName;

        m.plotdefaults.invisibleplot=default_plot_parameters.invisibleplot;
        m.plotdefaults.uicontrols=default_plot_parameters.uicontrols;
        m.pictures=default_plot_parameters.pictures;
    end
end

function m=savepng(m,NewDir)
    stagesuffix = makestagesuffixf( m.globalDynamicProps.currenttime );
    printfilename=fullfile(NewDir(1:end-4),'pngs',sprintf('E%s.png',stagesuffix));
    if exist(fullfile(NewDir(1:end-4),'pngs'))~=7
        mkdir(fullfile(NewDir(1:end-4),'pngs'));
        %             disp(sprintf('made directory %s',fullfile(NewDir(1:end-4),'pngs')));
    else
        %             disp(sprintf('directory %s exists',fullfile(NewDir(1:end-4),'pngs')));
    end
    if ispc
        m=leaf_plot(m,'invisibleplot',true,'uicontrols', false);
        disp('done leaf_plot');
        print(gcf,'-dpng',printfilename);
    else
        m=leaf_plot(m,'invisibleplot',true,'uicontrols', false);
        disp('done leaf_plot');
        print(gcf,'-dpng','-noui',printfilename);
    end
    disp(sprintf('written %s',printfilename))
end

function m=temporary_replacement_of_name(m,subdir,NewDir);
    temp_VideoName=m.globalProps.mov.VideoName(1:end-4);
    if exist(NewDir(1:end-4))~=7
        mkdir(NewDir(1:end-4));
        %         disp(sprintf('made directory %s',NewDir(1:end-4)));
    else
        %         disp(sprintf('directory %s exists',NewDir(1:end-4)));
    end
    m.globalProps.mov.VideoName=fullfile(temp_VideoName,subdir,'_000');
    tempdir=fullfile(NewDir(1:end-4),subdir);
    if exist(tempdir)~=7
        mkdir(tempdir);
        %         disp(sprintf('made directory %s',tempdir));
    else
        %         disp(sprintf('directory %s exists',tempdir));
    end
end

function [mov_VideoName,NewDir]=setup_movie_folder_and_names(Directory,ProjectName,mov_VideoName)
    if ~(exist(fullfile(Directory,ProjectName,'runs'))==7)
        mkdir(fullfile(Directory,ProjectName,'runs'));
    end
    if nargin<3
        mov_VideoName='';
    end
    disp(sprintf('setup_movie_folder_and_names:  mov_VideoName=%s',mov_VideoName));
    if isempty(mov_VideoName)
        DD=[];
        while isempty(DD) % generate unique directory name based on time
            moviename=[ProjectName,'-',datestr(now)];
            ind=findstr(moviename,' ');
            moviename(ind)='T';
            ind=findstr(moviename,':');
            moviename(ind)='-';
            rand('twister',floor(sum(1000*clock)));
            moviename=[moviename,'-R',num2str(floor(10*rand)),num2str(floor(10*rand)),num2str(floor(10*rand)),'R'];
            moviename=[moviename,'_000'];
            DD=dir(fullfile(Directory,ProjectName,'runs',moviename(1:end-4)));
            if isempty(DD)
                %             m.globalProps.mov.VideoName=moviename;
                mov_VideoName=moviename;
                NewDir=fullfile(Directory,ProjectName,'runs',mov_VideoName);
                mkdir(NewDir(1:end-4));
                DD=dir(fullfile(Directory,ProjectName,'runs',moviename(1:end-4)));
            end
        end
    else
        NewDir=fullfile(Directory,ProjectName,'runs',mov_VideoName);
        mkdir(NewDir(1:end-4));
    end
end

function keep_record_of_experiment(ProjectName,Directory,NewDir,inputparameters)
    cells=inputparameters;
    % save a copy of the command line that invoked this function
    if ischar(inputparameters)
        cmdstr = inputparameters;
        cells = {};
    else
        cmdstr='';
    end
    for i=1:2:length(cells)
        cmdstr=[cmdstr,'''',cells{i},'''',','];
        opt=cells{i+1};
        if iscell(opt)
            cmdstr=[cmdstr,'{'];
            for k=1:length(opt)
                if ischar(opt{k})
                    ch=['''',opt{k},''''];
                    cmdstr=[cmdstr,ch];
                else
                    ch=opt{k};
                    cmdstr=[cmdstr,num2str(ch),', '];
                end
            end
            cmdstr=[cmdstr,'},'];
        else
            if ischar(opt)
                ch=['''',opt,''''];
            else
                ch=num2str(opt);
            end
            cmdstr=[cmdstr,' ',num2str(ch),', '];
        end
    end
    h=fopen(fullfile(NewDir(1:end-4),'CommandLine.txt'),'w');
    fprintf(h,'%s',cmdstr);
    fclose(h);
    % save a copy of the interaction function with suffix .txt
    %remember hyphens are replaced by underscores
    ind=strfind(ProjectName,'-');
    ProjectNameU=lower(ProjectName);
    ProjectNameU(ind)='_';
    filefrom=fullfile(Directory,ProjectName,[ProjectNameU,'.m']);
    if exist(filefrom)==2
        copyfile(filefrom,...
            fullfile(NewDir(1:end-4),[ProjectNameU,'.txt']));
    end
end
