
% get current directry
cwd=pwd;

ROI_file=spm_get(1, '.mat', 'Select ROI.mat file to plot');

load(ROI_file);

time=[]; type=[]; vec_types=[]; pst=[];

if exist('plot_info.mat') == 2
  load plot_info;
end

%--------------------------------------
%  Collect details needed for averaging
%--------------------------------------

     %---- define the peristimulus time ----
      if exist('time')==0;

       num_scans=input('How many scans per trial: ');
       num_trials=input('How many trials in the experiment: ');
       pst=[1:num_scans];
      
       time=[];
        for i=1:num_trials     
         time=[time pst];  
        end

        %-Quick check
	fprintf('\nThere should have been ');
        fprintf('%8d',size(time,2));
        fprintf(' scans altogether in the experiment.\n')
        fprintf('\n');
      
      end 

      %-------- define condition types -------
      if exist('trial_type')==0;
     
         %-ask for names of the different trial types
         num_types=input('How many trial types (conditions): ');

             for i=1:num_types
              fprintf('Type'); fprintf('%2d',i);
              trial_type(i,:)=input([' - please enter name: '],'s');
             end
      end


      %-------- define trial types -----------
      if exist('vec_types')==0;
         
         fprintf('\n\nPlease, enter (the name of a matlab variable ');
         fprintf('representing) \na vector of length ');
         fprintf('%4d',num_trials);
         fprintf(' (or the number of trials). \n\nThe vector is needed in ');
         fprintf('order to define each trial as to the \ncondition ');
         fprintf('it belongs to. \nEnter 1 for '); fprintf(trial_type(1,:)); 
         fprintf(', 2 for '); fprintf(trial_type(2,:)); 
         fprintf(' and so forth.\n\n');
   
         vec_types=input('Vector defining the types: ');

      end

      %---------------     


      % replicate the vec_types to form an index vector, identifying
      % each scan as to which scondition it belongs to
      %-------------------------------
      type_idx = reshape(repmat(vec_types,num_scans,1),1,size(time,2));

      
      % Save plotting information in plot_info.mat
      cd(cwd);      
      save plot_info num* pst time trial_type vec_types type_idx;
        

  %=============================================
  %  start cycling through subjects if option=1      
  %=============================================
  for subject=1:avgROI.num_subj 
       
      avgROI.Y        = ROIdat{subject}.Y_mean;
      avgROI.Yc       = ROIdat{subject}.Yc_mean;
      avgROI.y        = ROIdat{subject}.y_mean;     
      avgROI.Yng      = ROIdat{1}.Yng_mean;

      %-----------------------------------------------------
      %-average Yc and y across time
      %-----------------------------------------------------

      for scan=1:num_scans                    
         for type=1:num_types

            t_avgROI.Y(scan,type)=...
                    mean(avgROI.Y  (find((time==scan)&(type_idx==type)))); 

            t_avgROI.Yc(scan,type)=...
                    mean(avgROI.Yc (find((time==scan)&(type_idx==type)))); 

            t_avgROI.y(scan,type)=...
                    mean(avgROI.y  (find((time==scan)&(type_idx==type)))); 

            t_avgROI.Yng(scan,type)=...
                    mean(avgROI.Yng(find((time==scan)&(type_idx==type)))); 
            
            
            % SEM(y) = StDev(y)/sqrt(n)
            %--------------------------

         end
      end

      %----------------------------------------------------- 
      %-plot fitted response, and adjusted data +/- SEM
      %----------------------------------------------------- 

      COL   = ['r','b','g','c','y','m'];

     figure
     hold on
     axis([0 max(pst) floor(min(min(t_avgROI.y)))  ceil(max(max(t_avgROI.y)))])

      for i=1:num_types
            plot(pst,t_avgROI.y(:,i),[COL(i) '-'],'LineWidth',3)
            legend(trial_type)
      end

      for i=1:num_types
            plot(pst,t_avgROI.Yc(:,i),[COL(i) '--'],'LineWidth',1)
      end

      title(['Fit & Adj data for ' ROIdat{subject}.name ','...
             ' for ' ROIdat{subject}.subj_path]);
%------------------------

     figure
     hold on
     axis([0 max(pst) floor(min(min(t_avgROI.Y)))  ceil(max(max(t_avgROI.Y)))])
      
      for i=1:num_types
       plot(pst,t_avgROI.Y(:,i),[COL(i) '-'],'LineWidth',2)
       legend(trial_type)
      end

      title(['Raw (scaled) data for ' ROIdat{1}.name]);
%--------------------------

    figure
    hold on
    axis([0 max(pst) floor(min(min(t_avgROI.Yng)))...
                      ceil(max(max(t_avgROI.Yng)))])

      for i=1:num_types
       plot(pst,t_avgROI.Yng(:,i),[COL(i) '-'],'LineWidth',3)
       legend(trial_type)
      end

      title(['Raw data for ' ROIdat{subject}.name ','...
             ' for ' ROIdat{subject}.subj_path]);
%---------------------------


   %-Create a structure with the trial-averaged timecourse for each subject
   % and for each condition 



  %================================================
  % end cycling through subjects, if this option
  %================================================
  end

