function AVL_batch_model4_spm5


%__________________________________________________________________________
% Batch script for AVL experiment:
%   - preprocessing
%       - realignment
%       - unwarp (optional)
%       - slice-timing (optional)
%       - spatial normalisation
%       - smoothing => 2 separate smoothing steps (4mm & 8mm)
%   - GLM estimation (event-related model)
%   - contrast definition & evaluation
% for a series of multi-session analyses.
% 3 basis functions and parametric modulation by time
%
% Klaas Enno Stephan & Hanneke den Ouden
% last changes: 06-12-2005
%
% If you want to modify this script for your own purposes, simply search for 
% and adapt all sections that have "#mod#" in their preceding comment.
%
% NB check beforehand whether onsets and RW are in the correct analysis
% folder!
%__________________________________________________________________________
%
%


% call SPM to have the graphics windows in place (#mod#)
%===========================================================================
addpath C:\Documents and Settings\mjgrol\MyDocuments\MATLAB\mei_spm5\spm5
addpath C:\Documents and Settings\mjgrol\MyDocuments\MATLAB\mei_spm5\spm_updates
spm fmri
%spm_defaults_AVL
global defaults
fs = filesep; % platform-specific file separator


% general variables (#mod#)
%===========================================================================
n_sess          = 1;
n_scans         = [360]; % scans per session
n_slices        = 34;
TR              = 3.22;
% for slice-timing
slice_order     = [1:1:n_slices];
ref_slice       = 17;
TR_slice        = 0.065;
TA              = (n_slices - 1) * TR_slice;
% SPM: time binning per TR
n_bins          = 16; % bins per TR
ref_bin         = 8; % reference bin
% names & paths
scans_prfx      = 'f';  
dir_base        = 'C:\Documents and Settings\mjgrol\My Documents\MATLAB\oefensets\fMRI_epoch_att2mot_ss\attention';
dir_analysis    = 'model';
dir_functional  = 'functional'; % base directory of functional scans (Analyze)
dir_protocol    = 'C:\Documents and Settings\mjgrol\My Documents\MATLAB\oefensets\fMRI_epoch_att2mot_ss\attention\protocols'; % dir where protocol files are written
sess_prfx       = ''; % prefix of session directories 'sess
protocol_prfx   = 'BatchProtocol_'; % prefix of protocol files
% preprocessing strategy
unwarp          = 1; % use unwarp?  1=yes, 0=no
rm_unwarped     = 0; % remove unwarped scans to save disk space? 1=yes, 0=no
st              = 0; % use slice-timing? 1=yes, 0=no
rm_st           = 0; % remove slice-timed scans to save disk space? 1=yes, 0=no
rp              = 0; % include realignment parameters in design matrix? 1=yes, 0=no


% subject-specific variables (#mod#)
%===========================================================================
name_subj       = {'M00587'};
% subject-specific start of analysis
start_analysis  = [4];
% 1 = realignment & unwarp
% 2 = slice-timing
% 3 = normalization & smoothing
% 4 = design matrix definition & parameter estimation
% 5 = define and compute contrasts
% 6 = smooth
% 7 = omit (subject has been analyzed)
% safety check
if length(name_subj) ~= length(start_analysis)
    disp('Error in variable definitions: different numbers of names and starting points!')
    return
end


% definition of experimental conditions (i.e. regressors in design matrix) 
% (#mod#)
%===========================================================================
cond_names      = {'Photic', 'Motion', 'Attention'};
cond_n          = length(cond_names);
cond_emptyRegr  = 0; % no. of conditions where empty regressors (zero events) can potentially occur




% loop through all subjects
%===========================================================================
k       = 1;


while (k<=length(name_subj)),
    
    err = 0;
    msg = {};
    
    % Set up some subject-specific things 
    % (unless this subject is already analyzed)
    % -------------------------------------------------------------------
    if start_analysis(k)<6
        
        % initialize subject-specific protocol file
        % -----------------------------------------
        cd (dir_protocol);
        protocol        = [protocol_prfx name_subj{k}];
        [fid,m]         = fopen([protocol '.txt'], 'wt+');
        
        % define subject-specific scan directories
        % ----------------------------------------
        for sess    = 1:n_sess,
            dir_scans{sess} = [dir_base fs name_subj{k} fs dir_functional];
        end
        
        % load subject-specific onsets (#mod#)
        % ------------------------------------
        try
            % create cell array cond_ons in which all conditions are stored
            cond_ons    = cell(n_sess,cond_n);
            % load onsets from external file
            cd ([dir_base fs name_subj{k} fs dir_analysis])
            load (['factors_' name_subj{k}]);
            % loop through sessions and assign onsets and durations to each condition
            for sess    = 1:n_sess,
                % onsets
                cond_ons{sess,1}  = phot;
                cond_ons{sess,2}  = mot;
                cond_ons{sess,3}  = att;

                % durations
                cond_dur{sess,1}  = 10;                      
                cond_dur{sess,2}  = 10;
                cond_dur{sess,3}  = 10;
              
            end
        catch
            err = 1;
            msg = logProc(msg,['Error loading condition vectors of subject ' name_subj{k}], fid);
        end
        
    end % start_analysis(k)<6 ?
    
    
    % realign & unwarp
    %===========================================================================
    if ~err & (start_analysis(k)<=1)
        
        % load Analyze files
        % NB:   paths (P) are stored session-wise, whereas V contains all images
        %       across all sessions
        %--------------------------------------------------------------------
        if ~err
            try
                % switch to scan directory of first session (to ensure that any spm2.ps file is written there)
                cd (dir_scans{1})
                % define cell array to store file and path names
                P   = cell(1,n_sess);
                for sess    = 1:n_sess,
                    % load scan file names into cell array 
                    % (NB: P{sess}(i,:) is the name of the i-th image)
                    P{sess} = spm_select('Files', dir_scans{sess}, [scans_prfx '*.img']);
                    disp(['Nbr of scans in session ' num2str(sess) ': ' num2str(size(P{sess},1))]);
                end
                V   = spm_vol(P);            
            catch
                err = 1;
                msg = logProc(msg,['Error loading Nifti files of subject ' name_subj{k}], fid);
            end
        end
        
        % Realignment (coregister only)
        %------------------------------------
        FlagsC  = struct(   'quality',1,...
            'fwhm',5,...
            'rtm',1,...
            'wrap',[0 1 0],...
            'interp',7);
        msg = logProc(msg,['Starting realignment (coregister only) for ' name_subj{k}], fid);
        try
            spm_realign(V,FlagsC);
        catch
            err = 1;
            msg = logProc(msg,['Error during realignment for subject ' name_subj{k}], fid);
        end
        
        % Unwarp (and reslice)
        %------------------------------------
        if unwarp
            
            msg = logProc(msg,['Starting unwarping (incl. reslicing) for ' name_subj{k}], fid);
            
            % parameters for estimation of unwarp parameters
            % NB: the current values correspond to the default parameters in spm_uw_estimate
            % NB: uw_est_flags.M points to the first image of the first session
            uw_est_flags = struct(  'order',           [12 12],...
                'sfP',             [],...
                'M',               V{1}(1,1).mat,...
                'regorder',        1,...
                'lambda',          1e5,...
                'jm',              0,...
                'fot',             [4 5],...
                'sot',             [],...
                'fwhm',            4,...
                'rem',             1,...
                'exp_round',       'Average',...
                'noi',             5,...
                'hold',            [1 1 1 0 1 0]);
            
            % Estimate unwarping parameters
            if ~err
                try
                    for sess    = 1:n_sess,
                        ds                  = spm_uw_estimate(P{sess},uw_est_flags);    % not sure what 'ds' stands for, just keeping Jesper's nomenclature...
                        uw_params(sess)     = ds;
                        [path,name,ext,ver] = fileparts(P{sess}(1,:));
                        pefile              = fullfile(path,[name '_uw.mat']);
                        save(pefile,'ds');
                    end
                catch
                    err = 1;
                    msg = logProc(msg,['Error estimating unwarp parameters for subject ' name_subj{k}], fid);
                end
            end
            
            % parameters for writing the unwarped images
            % NB: these correspond to the default parameters in spm_uw_apply,
            %     except for interp (here set to 7th order B-Spline)
            uw_write_flags = struct('mask',       1,...
                'mean',       1,...
                'interp',     7,...
                'wrap',       [0 1 0],...
                'which',      1,...
                'udc',        1);
            
            % Write unwarped images
            if ~err
                try
                    spm_uw_apply(uw_params,uw_write_flags);
                    %NB: no need to loop through sessions since spm_uw_apply
                    %    recognizes the number of sessions from the length of uw_params
                catch
                    err = 1;
                    msg = logProc(msg,['Error writing unwarped images of subject ' name_subj{k}], fid);
                end
            end
            
        end % unwarp?
        
    end % if-block for realign & unwarp section
    
    
    
    % slice-timing
    %===========================================================================
    if ~err & (start_analysis(k)<=2) & st
        
        msg = logProc(msg,['Starting slice-timing for ' name_subj{k}], fid);
        
        if unwarp   
            % load unwarped files
            %-----------------------------
            try
                P = cell(1,n_sess);
                for sess    = 1:n_sess,
                    P{sess} = spm_select('Files', dir_scans{sess}, ['u' scans_prfx '*.img']);
                end
            catch
                err = 1;
                msg = logProc(msg,['Error loading unwarped files of subject ' name_subj{k}], fid);
            end
        else
            % load realigned files
            %-----------------------------
            try
                P = cell(1,n_sess);
                for sess    = 1:n_sess,
                    P{sess} = spm_select('Files', dir_scans{sess}, [scans_prfx '*.img']);
                end
            catch
                err = 1;
                msg = logProc(msg,['Error loading realigned files of subject ' name_subj{k}], fid);
            end
        end
        % call slice-timing routine
        %--------------------------
        spm_slice_timing(P, slice_order, ref_slice, [TR_slice (TR-TA)]);
        
    end % slice-timing
    
    
    
    % Normalization & smoothing
    %=======================================
    
    if ~err & (start_analysis(k)<=3)
        
        if unwarp   
            % load unwarped files
            %-----------------------------
            try
                P = cell(1,n_sess);
                for sess    = 1:n_sess,
                    if st
                        P{sess} = spm_select('Files', dir_scans{sess}, ['au' scans_prfx '*.img']);
                    else
                        P{sess} = spm_select('Files', dir_scans{sess}, ['u' scans_prfx '*.img']);
                    end
                end
                V   = spm_vol(P);
                V   = cat(1,V{:}); % concatenate scans across sessions for normalization
            catch
                err = 1;
                msg = logProc(msg,['Error loading unwarped files of subject ' name_subj{k}], fid);
            end
            % Memory map the mean image
            try
                mean_img    = spm_select('Files', dir_scans{1}, ['meanu' scans_prfx '*.img']);
                Vm          = spm_vol(mean_img);
            catch
                err = 1;
                msg = logProc(msg,['Error loading mean image of subject ' name_subj{k}], fid);                
            end
        else
            % load realigned files
            %-----------------------------
            try
                P = cell(1,n_sess);
                for sess    = 1:n_sess,
                    if st
                        P{sess} = spm_select('Files', dir_scans{sess}, ['a' scans_prfx '*.img']);
                    else
                        P{sess} = spm_select('Files', dir_scans{sess}, [scans_prfx '*.img']);
                    end
                end
                V   = spm_vol(P);
                V   = cat(1,V{:}); % concatenate scans across sessions for normalization
            catch
                err = 1;
                msg = logProc(msg,['Error loading realigned files of subject ' name_subj{k}], fid);
            end
            % Memory map the first image of the time series
            try
                Vm  = spm_vol(P{1}(1,:));
            catch
                err = 1;
                msg = logProc(msg,['Error loading the first image of subject ' name_subj{k}], fid);                
            end                
        end
        
        
        % Estimate normalisation parameters
        %----------------------------------
        if ~err
            msg         = logProc(msg,['Starting normalisation for subject ' name_subj{k}], fid);
            matname     = [spm_str_manip(Vm.fname,'sd') '_sn.mat'];
            VG          = fullfile(spm('Dir'),'templates','EPI.nii');
            try
                % determine parameters
                params      = spm_normalise(VG,Vm,matname,'','',defaults.normalise.estimate);
                % mask
                msk         = spm_write_sn(V,params,defaults.normalise.write,'mask');
            catch
                err = 1;
                msg = logProc(msg,['Error estimating normalisation parameters (or masking) for subject ' name_subj{k}], fid);
            end
        end
        
        % Write normalised & smoothed (8 mm) images
        %------------------------------------------
        % NB: this part applies normalisation and smoothing in one step, using two different smoothing kernels.
        %     This is useful if one wants to analyze the data at both lower and higher spatial resolution.
        %        suw images  = 4 mm FWHM
        %        ssuw images = 8 mm FWHM
        if ~err
            % set voxel size and interpolation quality (#mod#)
            defaults.normalise.write.vox    = [3 3 3];
            defaults.normalise.write.interp = 7; % 7th degree B-spline interpolation
            % write normalised mean image
            msg = logProc(msg,['Writing normalised & smoothed images for subject ' name_subj{k}], fid); 
            try
                spm_write_sn(Vm,params,defaults.normalise.write,msk);
            catch
                err = 1;
                msg = logProc(msg,['Error writing normalised mean image for subject ' name_subj{k}], fid);
            end 
            % write normalised & smoothed images
            if ~err
                try
                    for i=1:length(V),
                        VO              = spm_write_sn(V(i),params,defaults.normalise.write,msk);
                        [pth,nam,ext]   = fileparts(V(i).fname);
                        fname           = fullfile(pth,['sw' nam ext]);
                        % smoothing - step 1 (suw images = 4 mm) (#mod#)
                        spm_smooth(VO,fname,4); 
                        % smoothing - step 2 (ssuw images = 8 mm) (#mod#)
                        fname           = fullfile(pth,['ssw' nam ext]);
                        spm_smooth(VO,fname,8);                 
                    end;
                    clear VO;
                catch
                    err = 1;
                    msg = logProc(msg,['Error writing normalised & smoothed images for subject ' name_subj{k}], fid);
                end
            end
        end
        
        % Remove unwarped & slice-timed images (which are no longer needed) if
        % desired by the user
        %-----------------------------------------------------------------
        if ~err & rm_unwarped
            msg = logProc(msg,['Deleting unwarped realigned images for subject ' name_subj{k}], fid);
            for sess    = 1:n_sess,
                cd (dir_scans{sess});
                delete (['u' scans_prfx '*.img']);
                delete (['u' scans_prfx '*.hdr']);
                delete (['u' scans_prfx '*.mat']);
            end
        end
        
        if ~err & rm_st
            msg = logProc(msg,['Deleting slice-timed images for subject ' name_subj{k}], fid);
            if unwarp
                prfx    = ['au' scans_prfx];
            else
                prfx    = ['a' scans_prfx];
            end
            for sess    = 1:n_sess,
                cd (dir_scans{sess});
                delete ([prfx '*.img']);
                delete ([prfx '*.hdr']);
                delete ([prfx '*.mat']);
            end
        end
        
    end    % if-block for normalisation & smoothing
    
    
    % define design matrix
    %===========================================================================
    
    % find empty regressors
    %---------------------------------------------------------------------------
    if start_analysis(k)<6
        emptyRegr = zeros (n_sess,1);
        for sess=1:n_sess,
            for c=1:cond_n
                if length(cond_ons{sess,c}) == 0
                    % empty regressor (zero events) encountered
                    emptyRegr(sess) = emptyRegr(sess) + 1;
                end
            end
        end
    end
    
    if ~err & (start_analysis(k)<=4)
        
        msg = logProc(msg,['Setting up design matrix for subject ' name_subj{k}], fid);
        
        % switch to analysis directory (to ensure all files are written there)
        try
            cd ([dir_base fs name_subj{k} fs dir_analysis]);
        catch
            err = 1;
            msg = logProc(msg,['Error switching to analysis directory of subject ' name_subj{k}], fid);
        end
        
        % number of scans and sessions & TR
        %---------------------------------------------------------------------------
        SPM.nscan          = ones(1,n_sess).*n_scans;
        SPM.xY.RT          = TR;
        
        % basis functions and timing parameters (#mod#)
        %---------------------------------------------------------------------------
        SPM.xBF.name       = 'hrf';
        SPM.xBF.order      = 1;                
        SPM.xBF.length     = 32;                % length in seconds
        SPM.xBF.T          = n_bins;            % number of time bins per scan
        SPM.xBF.T0         = ref_bin;               % middle slice/timebin
        SPM.xBF.UNITS      = 'scans';           % OPTIONS: 'scans'|'secs' for onsets
        SPM.xBF.Volterra   = 1;                 % OPTIONS: 1|2 = order of convolution
        
        % Trial specification: onsets, duration (in scans) and 
        % parameters for modulation  (#mod#)
        %---------------------------------------------------------------------------
        
        try
            % create cell array cond_ons in which all conditions that are parametrically 
            % modulated are stored (i.e. not conditions for target stimuli). 
            V    = cell(n_sess,(cond_n-2));
            dVdc = cell(n_sess,(cond_n-2));
            % load onsets from external file
            %cd ([dir_base fs name_subj{k} fs dir_analysis]);
            %load (['RW3_0_05_' name_subj{k}]);
            % loop through sessions and assign onsets and durations to each condition
        catch
            err = 1;
            msg = logProc(msg,['Error loading RW vectors of subject ' name_subj{k}], fid);
        end
        
        for sess=1:n_sess,
            cond = 1;
            
%             V{sess,1} = RW.V(sess).CSpAV;
%             V{sess,2} = RW.V(sess).CSpAv;
%             V{sess,3} = RW.V(sess).CSpaV;
%             V{sess,4} = RW.V(sess).CSpav;
%             V{sess,5} = RW.V(sess).CSmAV;
%             V{sess,6} = RW.V(sess).CSmAv;
%             V{sess,7} = RW.V(sess).CSmaV;
%             V{sess,8} = RW.V(sess).CSmav;
%             
%             dVdc{sess,1} = RW.dVdc(sess).CSpAV;
%             dVdc{sess,2} = RW.dVdc(sess).CSpAv;
%             dVdc{sess,3} = RW.dVdc(sess).CSpaV;
%             dVdc{sess,4} = RW.dVdc(sess).CSpav;
%             dVdc{sess,5} = RW.dVdc(sess).CSmAV;
%             dVdc{sess,6} = RW.dVdc(sess).CSmAv;
%             dVdc{sess,7} = RW.dVdc(sess).CSmaV;
%             dVdc{sess,8} = RW.dVdc(sess).CSmav;
            
            for c=1:cond_n
                if length(cond_ons{sess,c}) > 0
                    SPM.Sess(sess).U(cond).name      = {cond_names{c}};
                    SPM.Sess(sess).U(cond).ons       = cond_ons{sess,c};
                    SPM.Sess(sess).U(cond).dur       = cond_dur{sess,c};
                    % parametric modulation for first 8 conditions
%                     if c <= 8
%                         SPM.Sess(sess).U(cond).P(1).name = 'Param';
%                         SPM.Sess(sess).U(cond).P(1).P    = V{sess,c}';
%                         SPM.Sess(sess).U(cond).P(1).h    = 1; % 1st order polynomial expansion
%                         SPM.Sess(sess).U(cond).P(2).name = 'Param';
%                         SPM.Sess(sess).U(cond).P(2).P    = dVdc{sess,c}';
%                         SPM.Sess(sess).U(cond).P(2).h    = 1; % 1st order polynomial expansion
%                     else
                        % no parametric modulation for last two cond. (targets)
                        SPM.Sess(sess).U(cond).P.name = 'none';
%                         SPM.Sess(sess).U(cond).P(1).name = 'none';
%                         SPM.Sess(sess).U(cond).P(1).P = [];
%                         SPM.Sess(sess).U(cond).P(1).h = 0;
%                         SPM.Sess(sess).U(cond).P(2).name = 'none';                                            
%                         SPM.Sess(sess).U(cond).P(2).P = [];
%                         SPM.Sess(sess).U(cond).P(2).h = 0;                        
                        
%                     end
                end
                cond = cond+1;
            end
        end
        
        % User-specified regressors
        %--------------------------
        if rp
            % realignment parameters
            rnames = {'translX','translY','translZ','rotX','rotY','rotZ'};
            for sess = 1:n_sess
                fn = spm_select('Files',dir_scans{sess},'rp_*.txt');
                [r1,r2,r3,r4,r5,r6] = textread(fn,'%f%f%f%f%f%f');
                SPM.Sess(sess).C.C = [r1 r2 r3 r4 r5 r6];
                SPM.Sess(sess).C.name = rnames;
            end
        else
            % no user specified regressors 
            % (no need to include movement parameters when choosing to unwarp)
            for sess=1:n_sess
                SPM.Sess(sess).C.C = [];
                SPM.Sess(sess).C.name = {};
            end
        end
        
        % Define design matrix and save an additional copy (named
        % SPM_DesignMatrix.mat)
        %--------------------------------------------------------------------------
        try
            [SPM] = spm_fmri_design(SPM)
            save SPM_DesignMatrix SPM
        catch
            err = 1;
            msg = logProc(msg,['Error defining design matrix for subject ' name_subj{k}], fid);            
        end
        
        
        % Attach data to design matrix
        %===========================================================================
        
        if ~err
            msg = logProc(msg,['Attaching scans to design matrix for subject ' name_subj{k}], fid);
            
            % global normalization: OPTIONS:'Scaling'|'None'
            %---------------------------------------------------------------------------
            SPM.xGX.iGXcalc    = 'None';
            
            % low frequency confound: high-pass cutoff (secs) [Inf = no filtering]
            %---------------------------------------------------------------------------
            for sess=1:n_sess
                SPM.xX.K(sess).HParam = 128;
            end
            
            % intrinsic autocorrelations: OPTIONS: 'none'|'AR(1) + w'
            %-----------------------------------------------------------------------
            SPM.xVi.form       = 'AR(1)';
            
            % specify data by matrix of filenames (#mod#)
            %--------------------------------------------
            try
                tmpP        = cell(1,n_sess);
                for sess    = 1:n_sess,
                    % load scan file names into cell array (#mod#)
                    if unwarp
                        if st
                            tmpP{sess}  = spm_select('Files', dir_scans{sess}, (['snff' name_subj{k} '_*.img'])); %NB: change the prefix to 'swau' for manually preprocessed data
                        else
                            %filefilter = ['snff*.img'];
                            dirfiles = dir_scans{sess};
                            tmpP{sess} = spm_select('FPList', dirfiles, '^snff.*\.img$');
                            %fn = [repmat(dir_scans{sess},numel(fn(:,1)),1) strvcat(fn)];

                            %tmpP{sess}  = spm_select('Files', dir_scans{sess}, (['snff' name_subj{k} '_*.img'])); %NB: change the prefix to 'swu' for manually preprocessed data
%                             tmpP{sess}= spm_select('List',dir_scans{sess},['sswau' scans_prfx '*.img']);
                        end
                    else
                        if st
                            tmpP{sess}  = spm_select('Files', dir_scans{sess}, ['snff' name_subj{k} '_*.img']); %NB: change the prefix to 'swa' for manually preprocessed data
                        else
                            tmpP{sess}  = spm_select('Files', dir_scans{sess}, ['snff' name_subj{k} '_*.img']); %NB: change the prefix to 'sw' for manually preprocessed data
                        end
                    end
                end
                % concatenate scans across all sessions
                % (spm_fmri_spm_ui expects that all file names are in SPM.xY.P)
                SPM.xY.P    = cat(1,tmpP{:});
            catch
                err = 1;
                msg = logProc(msg,['Error loading preprocessed scans for subject ' name_subj{k}], fid);
            end
        end
        
        % Configure design matrix
        %-------------------------
        if ~err
            try
                SPM = spm_fmri_spm_ui(SPM);
            catch
                err = 1;
                msg = logProc(msg,['Error configuring design matrix for subject ' name_subj{k}], fid);
            end
        end
        
        
        % Estimate parameters
        %===========================================================================
        
        if ~err
            msg = logProc(msg,['Estimating parameters for subject ' name_subj{k}], fid);
            try
                SPM = spm_spm(SPM);
            catch
                err = 1;
                msg = logProc(msg,['Error estimating parameters for subject ' name_subj{k}], fid);
            end
        end
        
    end % if-block for design matrix definition & parameter estimation
    
    
    
    %     % Add contrasts
    %     %===========================================================================
    
    if ~err & (start_analysis(k)<=5)
        
        msg = logProc(msg,['Computing contrasts for subject ' name_subj{k}], fid);
        
        % switch to analysis directory and load SPM
        %------------------------------------------
        try
            clear SPM;
            cd ([dir_base fs name_subj{k} fs dir_analysis]);
            load SPM;            
        catch
            err = 1;
            msg = logProc(msg,['Error loading SPM.mat of subject ' name_subj{k}], fid);
        end
        
        % Create 1st contrast for 'effects of interest'
        % (necessary because the current version of spm_spm returns an
        % empty array for SPM.xCon)
        %------------------------------------------------------------------
        try
            iX0     = [SPM.xX.iB SPM.xX.iG];
        catch
            iX0     = [];
        end
        SPM.xCon        = spm_FcUtil('Set','effects of interest','F','iX0',iX0,SPM.xX.xKXs);
        
        % Add further contrasts
        % (#mod#)
        %------------------------------------------------------------------
        n_cntr      = length(SPM.xCon);                         % number of contrasts defined so far
        eb          = eye(SPM.xBF.order);                       % basis functions
        
        
        
        % Realignment parameters included in design matrix?
        if rp
            % define part of contrast vectors that represents movement parameters
            cntr_rp = zeros (1,6);
        else
            % contrast vectors do not need to take into account movement parameters
            cntr_rp = [];            
        end
        
        
        
        
        % NOTE: CONTRAST DEFINITIONS BELOW DO NOT ACCOUNT FOR EMPTY
        % REGRESSORS OR REALIGNMENT PARAMETERS !!!!
        
        c_name{1}   = 'F: main effect: motion';
        c_wgt{1}    = [0 1 0]; 
        c_type{1}   = 'F';
        
        c_name{2}   = 'T: main effect: attention';
        c_wgt{2}    = [0 0 1];
        c_type{2}   = 'T';
  
        c_name{3}   = 'F: main effect: photic';
        c_wgt{3}    = [1 0 0];
        c_type{3}   = 'F';
        
        % save contrast definitions in SPM.xCon                
        try
            for c = 1:length(c_name)
                cw                      = [c_wgt{c} zeros(size(c_wgt{c},1), n_sess)]';  % pad with zeros for constants & transpose
                SPM.xCon(c+n_cntr)      = spm_FcUtil('Set',c_name{c},c_type{c},'c',cw,SPM.xX.xKXs);
            end
        catch
            err = 1;
            msg = logProc(msg,['Error setting up contrasts for subject ' name_subj{k}], fid);
        end
        
        
        % Evaluate contrasts
        %---------------------------------------------------------------------------
        if ~err
            try
                spm_contrasts(SPM);
            catch
                err = 1;
                msg = logProc(msg,['Error evaluating contrasts for subject ' name_subj{k}], fid);
            end
        end
        
        if ~err
            msg = logProc(msg,['Analysis successfully completed for subject ' name_subj{k}], fid);
            msg = logProc(msg,['====================================================='], fid);
        end
        
    end     % if-block for contrast definition & evaluation
    
    
    % Smooth con images for RFX by 3-D Gaussian of 6mm FWHM. 'con' are the
    % conditions that have to be smoothed. 
    % =====================================================================
    if ~err & (start_analysis(k)<=6)
        for con = [3]
            strin = sprintf('%s\\con_%04d.img',fullfile(dir_base, name_subj{k}, dir_analysis),con);
            strout = sprintf('%s\\scon_%04d.img',fullfile(dir_base, name_subj{k}, dir_analysis),con);
            spm_smooth(strin,strout,[6 6 6]);
        end
    end
    
    % Switch to next subject
    %=======================================
    % save protocol file for current subject
    if start_analysis(k)<6
        fclose(fid);
    end
    % increase k
    k   = k + 1;
    % set err to zero and clear major variables
    err = 0;
    clear SPM V Vm P tmpP;
    
end     % end of main loop


% save protocol for all subjectS in .mat file
%===========================================


if ~isempty(find(start_analysis>6))
    cd (dir_protocol);
    eval(['save ' protocol ' msg;']);
end
return


%_____________________________________
% used internally


function [msg] = logProc (msg,proc,fid)


% display in command window
disp(proc);


% write into text file
fprintf(fid, [proc '\n']);


% store in msg cell array 
msg = {msg{:},proc};

return 
