function [tw, cmpdnames, time_vals, span_vals] = twui_ave(num_runs, first_runid, output_dir, parsed_partname, analysis, plot_info, ssflag, binsize)

% 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



%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

deltat = analysis.timeend - analysis.timestart;
if binsize == 0
    numbins = 1;
else
    numbins = ceil(deltat/binsize);
end

cur_runid = first_runid;

tw = [];
time = [];
time_vals = [];
span_vals = [];
temp_data_time = [];
temp_data_cmpds = [];

df = findstr(parsed_partname,'_data');

for i = 1:num_runs
   this_filename = [parsed_partname '_' num2str(cur_runid)];
%if work with parsed data
data_time = [];
data_cmpds = [];
if (isempty(df))
   load(fullfile(output_dir, this_filename), 'data_time', 'data_cmpds')
else  
    [data_time, data_cmpds] = read_snapshot(fullfile(output_dir, this_filename));
     snapshot_interval =  plot_info.snapshot_interval;
            if (binsize < snapshot_interval) && i == 1
%                error( ' Time binsize is smaller then snapshot interval!!! ')
                warning( ' Time binsize is smaller then snapshot interval!!! ')
            end
            
end
          
           
    for k = 1:length(c_index)
        
        cmpd_idx = c_index(k);
        t_index_lower = [];
        t_index_upper = [];
        outval_idx = 0;
        
        if ssflag
            temp_data_time = data_time;
            temp_data_cmpds = data_cmpds(:,cmpd_idx);

        else
            temp_data_time = data_time{cmpd_idx};
            temp_data_cmpds = data_cmpds{cmpd_idx};
        end

        d_time = analysis.timestart;
        
        % find the overall beginning and end indices for preallocation of
        % time, time_vals, and span_vals
        
        t_interval = analysis.timeend - d_time;

        for j = 2:length(temp_data_time)
            if (temp_data_time(j) - d_time > 0)
                t_index_lower = j-1;
                break
            end
        end    
        if isempty(t_index_lower)
            t_index_lower = length(temp_data_time);
        end

        for j = t_index_lower:length(temp_data_time)
            if (temp_data_time(j) - d_time >= t_interval)
                t_index_upper = j-1;
                break    
            end  
        end   
        if isempty(t_index_upper)
            t_index_upper = length(temp_data_time);
        end
        
        tot_num_ts = t_index_upper-t_index_lower+1;
        temp_time = zeros(tot_num_ts,1);
        temp_time_vals = zeros(tot_num_ts,1);
        temp_span_vals = zeros(tot_num_ts,1);;
        
        
        for bi = 1:numbins
            
            twa = 0;
            
            if bi == numbins
                t_interval = analysis.timeend - d_time;
            else
                t_interval = binsize;
            end
            
            if bi ~= 1
                t_index_lower = t_index_upper;
                t_index_upper = [];
            end
                
            for j = t_index_lower:length(temp_data_time)
                if (temp_data_time(j) - d_time >= t_interval)
                    t_index_upper = j-1;
                    break    
                end  
            end   
            if isempty(t_index_upper)
                t_index_upper = length(temp_data_time);
            end




            index = t_index_lower;          
            num_ts = t_index_upper-t_index_lower+1;

            if binsize == 0
                t_index_upper = t_index_lower;
            end
            
            if t_index_upper ~= t_index_lower
                temp_time(outval_idx + 1) = temp_data_time(t_index_lower+1) - d_time;
                temp_time_vals(outval_idx + 1) = d_time;
                temp_time(outval_idx + num_ts) = d_time + t_interval - temp_data_time(t_index_upper);
                temp_time_vals(outval_idx + num_ts) = temp_data_time(t_index_upper);

                for n = 2:(t_index_upper - t_index_lower)
                    temp_time(outval_idx + n) = temp_data_time(n+t_index_lower) - temp_data_time(n+t_index_lower-1);
                    temp_time_vals(outval_idx + n) = temp_data_time(t_index_lower+n-1);
                end

                for n = outval_idx + 1:outval_idx + num_ts
                    temp_span_vals(n) = temp_data_cmpds(index);
                    index = t_index_lower + n - outval_idx;
                    twa = twa + temp_time(n)/(t_interval)*temp_span_vals(n);
                end     

            else              
                temp_time_vals(outval_idx + 1) = d_time;
                temp_span_vals(outval_idx + 1) = temp_data_cmpds(t_index_lower);
                twa = temp_data_cmpds(t_index_lower);
            end           
            

            tw(i,k,bi) = twa;
            
            d_time = d_time + binsize;
            outval_idx = outval_idx + num_ts;


        end % numbins
        
%        time{i,k} = temp_time;
        time_vals{i,k} = temp_time_vals;
        span_vals{i,k} = temp_span_vals;

    end % length(c_index)
    cur_runid = cur_runid + 1;
end % num_runs
        

