function [tw, cmpdnames] = tw_ave(num_runs, first_runid, output_dir, parsed_partname, analysis, plot_info)

% computes the time weighted average from the the parsed output files specified by
% output_dir and parsed_partname, and for the compounds and times specified
% in analysis


%switch lower(analysis.type)
    
%    case 'average'
        

        %load('s_data.mat');
        %load('analysis_data.mat');
%        savefile = fullfile(output_dir, 'twa.mat');

d_time = analysis.timestart;
t_interval = analysis.timeend - d_time;

c_index = [];
if strcmp(analysis.cmpd{1}, 'all')
    c_index = [1:length(plot_info.monitor_cmpds)];
else
    for j = 1:length(analysis.cmpd)
        for i = 1:length(plot_info.monitor_cmpds)       
            if strcmp(analysis.cmpd{j}, plot_info.monitor_cmpds(i)) == 1
               c_index(j) = i;             
               break
           end
        end
    end
end

cmpdnames = [];
for k = 1:length(c_index)
    cmpdnames{k} = plot_info.monitor_cmpds{c_index(k)};
end

cur_runid = first_runid;

for i = 1:num_runs

    this_filename = [parsed_partname '_' num2str(cur_runid)];
    load(fullfile(output_dir, this_filename), 'data_time', 'data_cmpds')

    for k = 1:length(c_index)

        twa = [];
        time = [];
        span_vals = [];

        for j = length(data_time{c_index(k)}): -1:1            
%                if (data_time{c_index}(j) - d_time <= 0)
            if (data_time{c_index(k)}(j) - d_time < t_interval)
                t_index_upper = j;
                break    
            end  
        end   
        for j = t_index_upper:-1:1     
%                if (abs(data_time{c_index}(j) - d_time) >= t_interval)
            if (data_time{c_index(k)}(j) - d_time < 0)
                t_index_lower = j;
                break
            end
        end    

        index = t_index_lower;            

        if t_index_upper ~= t_index_lower
%                time{i}(1) = (d_time - t_interval) - data_time{c_index}(t_index_lower);            
            time{i}(1) = data_time{c_index(k)}(t_index_lower+1) - d_time;
%                time{i}(t_index_upper-t_index_lower) = d_time - data_time{c_index,i}(t_index_upper);
            time{i}(t_index_upper-t_index_lower+1) = d_time + t_interval - data_time{c_index(k)}(t_index_upper);

            for n = 1:(t_index_upper-t_index_lower+1)      
                span_vals{i}(n) = data_cmpds{c_index(k)}(index);
                index = t_index_lower + n;
            end     

            for n = 2:(t_index_upper - t_index_lower)
                time{i}(n) = data_time{c_index(k)}(n+t_index_lower) - data_time{c_index(k)}(n+t_index_lower-1);                                           
            end            
        else              
%                time{i}(1) = (d_time - t_interval) - data_time{c_index}(t_index_lower);        
            time{i}(1) = t_interval;        
%                span_vals{i}(1) = data_cmpds{c_index}(1);
            span_vals{i}(1) = data_cmpds{c_index(k)}(t_index_lower);
        end           

        for j = 1:length(span_vals{i})           
            twa{i}(j) = time{i}(j)/(sum(time{i}))*span_vals{i}(j);
        end

        tw{k}(i) = sum(twa{i});
%            sd(i) = std(span_vals{i});
%                save(savefile, 'tw');

    end
    cur_runid = cur_runid + 1;
end
        
%        mean_tw = zeros(1,num_runs/100);
%        mean_sd = zeros(1,num_runs/100);

%        for j = 0:num_runs/100-1        
%        for i = 1:num_runs
%           mean_tw = tw(i) + mean_tw;
%           mean_sd = sd(i) + mean_sd;
%        end

%        for k = 1:length(c_index)
%            mean_tw(k) = mean(tw{k})
%            std_tw(k) = std(tw{k})
%            figure
%            hist(tw{k})
%            hleg = legend(plot_info.monitor_cmpds{c_index(k)});
%            if this_matlab_version == 7;
%                set(hleg, 'Interpreter','none');
%            end
%        end
%        drawnow
%        end
        

%        for k = 0:num_runs/100-1
%        for j = 1:100
%            b{k+1}(j) = tw(j + k*100);
%        end          
%        hist(b{k+1})    
%        hist(tw)    
%        pause
%        end
        
        
%end
