function simres = grsim_simsweeppar(vlist, flist, olist, Nh, simfunc, parfunc)
%
% simres = grsim_simsweeppar(vlist, flist, olist, Nh, simfunc, parfunc)
%
% [Parallel version of grsim_simsweep.]
% Run the grading simulation function 'simfunc', to simulate
% each point n whose amplitude (v), frequency (f) and other
% (o [string]) parameters are given as vlist(n), flist(n)
% and olist{n}.
% Return time-domain data in the 'simres' structure, along
% with transformed (sine-based) frequency-domain data up to
% the Nh'th harmonic.
%
% Nathaniel, 2008-05.

% a primitive attempt at finding the directory containing
% this script, i.e. like `dirname $0`
[grsimd, myname, myext] = fileparts(mfilename('fullpath'));

% set a new directory for the temporary files for batching,
% outside the normal $HOME directory to avoid backup of junk
clk = clock;
udate = sprintf('%4d-%02d-%02d_%02d%02d%02d', ...
    clk(1), clk(2), clk(3), clk(4), clk(5), round(clk(6)) );
bdir = ['/home/public/nt_sim/grsim__', udate];

% make or copy necessary parts in the batch working-dir
unix(['mkdir -pv ', bdir, '/logs']);
unix(['cp ', grsimd, '/batchfiles/*', ' ', bdir, '/']);
unix(['mkdir -pv ', bdir, '/mfiles']);
unix(['cp ', grsimd, '/*.m', ' ', bdir, '/mfiles/']);
disp(' ')

% load default settings for simulation:
% geometry, material, solution-parameters
[geom, mtrl, solp] = parfunc();


% useful values from the input
N = numel(vlist);
Nt = solp.Nkeep*solp.Npph*Nh;

% allow a null olist to be made automatically by
% passing anything other than a cell array (e.g. 0, '')
if ~iscell(olist),
    olist = {''};
    for n=1:N,
        olist = {olist{:},''};
    end;
    olist = olist.';
end

% empty output arrays
realt = NaN*ones(N,1);
cput = NaN*ones(N,1);
tpoints = NaN*ones(N,Nt);
v_t = NaN*ones(N,Nt);
i_t = NaN*ones(N,Nt);
v_f = NaN*ones(N,Nh);
i_f = NaN*ones(N,Nh);


% set up for each parameter setting
disp('Making a directory and configuration data for each job:');
for n=1:N,

    % reload base settings for each
    [geom, mtrl, solp] = parfunc();

    v = vlist(n);
    f = flist(n);
    eval(olist{n});

    T = 1/f;
    dT = T/(solp.Npph*Nh);

    solp.vsup = sprintf('%5.3e*sin(2*pi*%5.3e*t)', v, f);
    fprintf('Point % 3d/%d:   %4.1fkV %1.2eHz\n', n, N, v/1e3, f);
    
    % a list of time-points to simulate for, with 'skipped'
    % (equilibrium-reaching) cycles at negative time; for a
    % later t>=0 check, it is made certain that there exists
    % a point exactly at zero
    solp.tsoln = [ (-solp.Nskip*T):dT:(-dT), 0:dT:(solp.Nkeep*T-dT) ];
    tidx = find(solp.tsoln>=0);
    
    % save whole settings lists for each point
    settings.geom{n} = geom;
    settings.mtrl{n} = mtrl;
    settings.solp{n} = solp;

    % create settings mat-file and job m-file, in job directory
    ndir = sprintf('%s/%06d', bdir, n);
    unix(['mkdir ', ndir]);
    save('-mat', [ndir,'/in.dat'], 'geom', 'mtrl', 'solp' );
    mscript = fopen([ndir,'/runme.m'], 'w');
    fprintf(mscript, [ ...
        '%% %06d\n', ...
        'maxNumCompThreads(1);\n', ...
        'load(''-mat'', ''in.dat'');\n', ...
        'path(''%s'',path);\n', ...
        '[tmpres, tmpfem] = %s(geom, mtrl, solp);\n', ...
        'save(''-mat'', ''out.dat'', ', ...
        '''tmpres'', ''geom'', ''mtrl'', ''solp'');\n', ...
        'exit\n' ],  n, [bdir,'/mfiles/'], simfunc);
    fclose(mscript);

    % (highly ugly): the 1dml simulation creates this
    % extra file when it is run, to be the pde function;
    % the file must exist before the 1dml function is
    % run, else it isn't found on the path; so, make a 
    % blank file
    unix(['echo "void" >', ndir, '/nlc_pde_tmp.m']);

end
disp(' ');

% Run batching scripts
disp('Starting jobs on all hosts and processors listed:');
unix([bdir,'/start_distributed_jobs']);
disp(' ');
% Wait for all jobs to finish
disp('Waiting for all jobs to finish...');
unix([bdir,'/wait_till_finished']);


% retrieve results for each parameter setting
disp(' ');
disp('Collecting output data from all job directories:');
clear simres;
for n=1:N,

    try

        ndir = sprintf('%s/%06d', bdir, n);
        clear tmpres geom mtrl solp;
        load('-mat', [ndir,'/out.dat']);

        cput(n)  = tmpres.tcpu;
        realt(n) = tmpres.treal;

        fprintf('%s:  % 3d/%d:  %4.1f s real, %4.1f s cpu\n', ...
            simfunc, n, N, realt(n), cput(n) );

        tidx = find(solp.tsoln>=0);

        % time-domain results into array of all
        tpoints(n,:) = solp.tsoln(tidx);
        v_t(n,:) = tmpres.v_t(tidx);
        i_t(n,:) = tmpres.i_t(tidx);
        
        % calculate frequency-domain results
        v_f(n,:) = grsim_fft(v_t(n,:), solp.Nkeep, Nh, 0);
        i_f(n,:) = grsim_fft(i_t(n,:), solp.Nkeep, Nh, 0);
        shift_v = -angle(v_f(n,1));
        v_f(n,:) = grsim_tshift(v_f(n,:), shift_v, 0);
        i_f(n,:) = grsim_tshift(i_f(n,:), shift_v, 0);

    catch

        disp('!! PROBLEM with retrieving simulation data');
        disp(['iteration ', num2str(n), ' in ', mfilename]);
        if exist('lasterr'); disp(lasterr); end

    end

end
disp(' ');

% output as a structure
simres.simfunc = simfunc;
simres.Nh = Nh;
simres.vlist = vlist;
simres.flist = flist;
simres.olist = olist;
simres.settings = settings;
simres.cput = cput;
simres.realt = realt;
simres.t = tpoints;
simres.v_t = v_t;
simres.i_t = i_t;
simres.v_f = v_f;
simres.i_f = i_f;