function write_fsf(fsf)
% write_fsf(fsf);
%
% Writes a FEAT structure file (.fsf) for use with FSL's FEAT analysis package
%
% fsf is structure with all of the parameters necessary to write out the .fsf
% file. The default field values are set in create_fsf.m
%
% NOTE: This is still a work in progress. 
% Fortunately, the FEAT GUI will read a partially completed setup file.  
% In effect, it should be possible to create a template for any given analysis
% and then fill a few fields with pointers to the subject specific files.

% 04/22/05 Petr Janata
% 09/05/05 PJ - completed first functional version

% Open the .fsf file based on information in the structure

fid = fopen(fullfile(fsf.fsldir,'design.fsf'),'wt');

fsf_str = sprintf([ ...
      '\n# FEAT version number' ...
      '\nset fmri(version) %1.2f\n'], fsf.version);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([ ...
      '\n# Are we in MELODIC?' ...
      '\nset fmri(inmelodic) 0\n']);
fprintf(fid,'%s',fsf_str);

fsf_str = sprintf([...
      '\n# Analysis level' ...
      '\n# 1 : First-level analysis' ...
      '\n# 2 : Higher-level analysis' ...
      '\nset fmri(level) %d\n'], fsf.level);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Which stages to run' ...
      '\n# 0 : No first-level analysis (registration and/or group stats only)' ...
      '\n# 7 : Full first-level analysis' ...
      '\n# 1 : Pre-Stats' ...
      '\n# 3 : Pre-Stats + Stats' ...
      '\n# 2 :             Stats' ...
      '\n# 6 :             Stats + Contrasts, Thresholding, Rendering' ...
      '\n# 4 :                     Contrasts, Thresholding, Rendering' ...
      '\nset fmri(analysis) %d\n'], fsf.analysis);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.4
  fsf_str = sprintf([...
        '\n# Delay before starting (hours)' ...
        '\nset fmri(delay) %d\n'], fsf.delay);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Use relative filenames' ...
      '\nset fmri(relative_yn) %d\n'], fsf.relative_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Balloon help' ...
      '\nset fmri(help_yn) %d\n'], fsf.help_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Run Featwatcher' ...
      '\nset fmri(featwatcher_yn) %d\n'], fsf.featwatcher_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Cleanup first-level standard-space images' ...
      '\nset fmri(sscleanup_yn) %d\n'], fsf.sscleanup_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Output directory' ...
      '\nset fmri(outputdir) "%s"\n'], fsf.outputdir);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# TR(s)' ...
      '\nset fmri(tr) %1.2f\n'], fsf.tr);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Total volumes' ...
      '\nset fmri(npts) %d\n'], fsf.npts);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Delete volumes' ...
      '\nset fmri(ndelete) %d\n'], fsf.ndelete);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf([...
      '\n# Perfusion tag/control order' ...
      '\nset fmri(tagfirst) %d\n'], fsf.tagfirst);
  fprintf(fid,'%s', fsf_str);
end
  
fsf_str = sprintf([...
      '\n# Number of first-level analyses' ...
      '\nset fmri(multiple) %d\n'], fsf.multiple);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Higher-level input type' ...
      '\n# 1 : Inputs are lower-level FEAT directories' ...
      '\n# 2 : Inputs are cope images from FEAT directories' ...
      '\nset fmri(inputtype) %d\n'], fsf.inputtype);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Carry out pre-stats processing?' ...
      '\nset fmri(filtering_yn) %d\n'], fsf.filtering_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Brain/background threshold, %%' ...
      '\nset fmri(brain_thresh) %d\n'], fsf.brain_thresh);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf([...
      '\n# Critical z for design efficiency calculation' ...
      '\nset fmri(critical_z) %1.2f\n'], fsf.critical_z);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
      '\n# Noise level' ...
      '\nset fmri(noise) %1.6f\n'], fsf.noise);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
      '\n# Noise AR(1)' ...
      '\nset fmri(noisear) %1.6f\n'], fsf.noisear);
  fprintf(fid,'%s', fsf_str);
end
  
fsf_str = sprintf([...
      '\n# Post-stats-only directory copying' ...
      '\n# 0 : Overwrite original post-stats results' ...
      '\n# 1 : Copy original FEAT directory for new Contrasts, Thresholding, Rendering' ...
      '\nset fmri(newdir_yn) %d\n'], fsf.newdir_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Slice timing correction' ...
      '\n# 0 : None' ...
      '\n# 1 : Regular up (0, 1, 2, 3, ...)' ...
      '\n# 2 : Regular down' ...
      '\n# 3 : Use slice order file' ...
      '\n# 4 : Use slice timings file' ...
      '\n# 5 : Interleaved (0, 2, 4 ... 1, 3, 5 ... )' ...
      '\nset fmri(st) %d\n'], fsf.st);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Slice timings file' ...
      '\nset fmri(st_file) "%s"\n'], fsf.st_file);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Motion correction' ...
      '\n# 0 : None' ...
      '\n# 1 : MCFLIRT' ...
      '\nset fmri(mc) %d\n'], fsf.mc);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Spin-history (currently obsolete)' ...
      '\nset fmri(sh_yn) %d\n'], fsf.sh_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# BET brain extraction' ...
      '\nset fmri(bet_yn) %d\n'], fsf.bet_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Spatial smoothing FWHM (mm)' ...
      '\nset fmri(smooth) %d\n'], fsf.smooth);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Intensity normalization' ...
      '\nset fmri(norm_yn) %d\n'], fsf.norm_yn);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf([...
      '\n# Perfusion subtraction' ...
      '\nset fmri(perfsub_yn) %d\n'], fsf.perfsub_yn);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Highpass temporal filtering' ...
      '\nset fmri(temphp_yn) %d\n'], fsf.temphp_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Lowpass temporal filtering' ...
      '\nset fmri(templp_yn) %d\n'], fsf.templp_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# MELODIC ICA data exploration' ...
      '\nset fmri(melodic_yn) %d\n'], fsf.melodic_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Carry out main stats?' ...
      '\nset fmri(stats_yn) %d\n'], fsf.stats_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Carry out prewhitening?' ...
      '\nset fmri(prewhiten_yn) %d\n'], fsf.prewhiten_yn);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf([...
      '\n# Add motion parameters to model' ...
      '\n# 0 : No' ...
      '\n# 1 : Yes' ...
      '\nset fmri(motionevs) %d\n'], fsf.motionevs);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
      '\n# Robust outlier detection in FLAME?' ...
      '\nset fmri(robust_yn) %d\n'], fsf.robust_yn);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Higher-level modelling' ...
      '\n# 3 : Fixed effects' ...
      '\n# 0 : Mixed Effects: Simple OLS' ...
      '\n# 2 : Mixed Effects: FLAME (stage 1 only)' ...
      '\n# 1 : Mixed Effects: FLAME (full)' ...
      '\nset fmri(mixed_yn) %d\n'], fsf.mixed_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Number of EVs' ...
      '\nset fmri(evs_orig) %d' ...
      '\nset fmri(evs_real) %d\n'], fsf.evs_orig,fsf.evs_real);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf('set fmri(evs_vox) %d\n', fsf.evs_vox);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Number of contrasts' ...
      '\nset fmri(ncon_orig) %d' ...
      '\nset fmri(ncon_real) %d\n'], fsf.ncon_orig, fsf.ncon_real);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Number of F-tests' ...
      '\nset fmri(nftests_orig) %d' ...
      '\nset fmri(nftests_real) %d\n'], fsf.nftests_orig, fsf.nftests_real);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Add constant column to design matrix? (obsolete)' ...
      '\nset fmri(constcol) %d\n'], fsf.constcol);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Carry out post-stats steps?' ...
      '\nset fmri(poststats_yn) %d\n'], fsf.poststats_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Pre-threshold masking?' ...
      '\nset fmri(threshmask) "%s"\n'], fsf.threshmask);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Thresholding' ...
      '\n# 0 : None' ...
      '\n# 1 : Uncorrected' ...
      '\n# 2 : Voxel' ...
      '\n# 3 : Cluster' ...
      '\nset fmri(thresh) %d\n'], fsf.thresh);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# P threshold' ...
      '\nset fmri(prob_thresh) %1.4f\n'], fsf.prob_thresh);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Z threshold' ...
      '\nset fmri(z_thresh) %1.2f\n'], fsf.z_thresh);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Z min/max for colour rendering' ...
      '\n# 0 : Use actual Z min/max' ...
      '\n# 1 : Use preset Z min/max' ...
      '\nset fmri(zdisplay) %d\n'], fsf.zdisplay);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Z min in colour rendering' ...
      '\nset fmri(zmin) %d\n'], fsf.zmin);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Z max in colour rendering' ...
      '\nset fmri(zmax) %d\n'], fsf.zmax);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Colour rendering type' ...
      '\n# 0 : Solid blobs' ...
      '\n# 1 : Transparent blobs' ...
      '\nset fmri(rendertype) %d\n'], fsf.rendertype);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Background image for higher-level stats overlays' ...
      '\n# 1 : Mean highres' ...
      '\n# 2 : First highres' ...
      '\n# 3 : Mean functional' ...
      '\n# 4 : First functional' ...
      '\n# 5 : Standard space template' ...
      '\nset fmri(bgimage) %d\n'], fsf.bgimage);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf([...
      '\n# Create time series plots' ...
      '\nset fmri(tsplot_yn) %d\n'], fsf.tsplot_yn);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Registration?' ...
      '\nset fmri(reg_yn) %d\n'], fsf.reg_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([... 
      '\n# B0 fieldmap unwarping?' ...
      '\nset fmri(regunwarp_yn) %d\n'], fsf.regunwarp_yn);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf([...
      '\n# EPI dwell time (ms)' ...
      '\nset fmri(dwell) %d\n'], fsf.dwell);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
      '\n# EPI TE (ms)' ...
      '\nset fmri(te) %d\n'], fsf.te);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
      '\n# % Signal loss threshold' ...
      '\nset fmri(signallossthresh) %d\n'], fsf.signallossthresh);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
      '\n# Unwarp direction' ...
      '\nset fmri(unwarp_dir) %s\n'], fsf.unwarp_dir);
  fprintf(fid,'%s', fsf_str);
end

if fsf.version == 5.4
  fsf_str = sprintf([...
        '\n# Dwell/Asymmetry ratio ' ...
        '\nset fmri(dwellasym) %d\n'], fsf.dwellasym);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Registration to initial structural' ...
      '\nset fmri(reginitial_highres_yn) %d\n'], fsf.reginitial_highres_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Search space for registration to initial structural' ...
      '\n# 0   : No search' ...
      '\n# 90  : Normal search' ...
      '\n# 180 : Full search' ...
      '\nset fmri(reginitial_highres_search) %d\n'], fsf.reginitial_highres_search);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Degrees of Freedom for registration to initial structural' ...
      '\nset fmri(reginitial_highres_dof) %d\n'], fsf.reginitial_highres_dof);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.4
  fsf_str = sprintf([...
        '\n# Do nonlinear registration to initial structural?' ...
        '\nset fmri(reginitial_highres_nonlinear_yn) %d\n'], fsf.reginitial_highres_nonlinear_yn);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Registration to main structural' ...
      '\nset fmri(reghighres_yn) %d\n'], fsf.reghighres_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Search space for registration to main structural' ...
      '\n# 0   : No search' ...
      '\n# 90  : Normal search' ...
      '\n# 180 : Full search' ...
      '\nset fmri(reghighres_search) %d\n'], fsf.reghighres_search);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Degrees of Freedom for registration to main structural' ...
      '\nset fmri(reghighres_dof) %d\n'], fsf.reghighres_dof);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.4
  fsf_str = sprintf([...
        '\n# Do nonlinear registration to main structural?' ...
        '\nset fmri(reghighres_nonlinear_yn) %d\n'], fsf.reghighres_nonlinear_yn);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Registration to standard image?' ...
      '\nset fmri(regstandard_yn) %d\n'], fsf.regstandard_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Standard image' ...
      '\nset fmri(regstandard) "%s"\n'], fsf.regstandard);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Search space for registration to standard space' ...
      '\n# 0   : No search' ...
      '\n# 90  : Normal search' ...
      '\n# 180 : Full search' ...
      '\nset fmri(regstandard_search) %d\n'], fsf.regstandard_search);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  fsf_str = sprintf([...
      '\n# Control nonlinear warp field resolution' ...
      '\nset fmri(regstandard_nonlinear_warpres) %d\n'],...
      fsf.regstandard_nonlinear_warpres);
  fprintf(fid,'%s', fsf_str);
end

fsf_str = sprintf([...
      '\n# Degrees of Freedom for registration to standard space' ...
      '\nset fmri(regstandard_dof) %d\n'], fsf.regstandard_dof);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Do nonlinear registration to standard space?' ...
      '\nset fmri(regstandard_nonlinear_yn) %d\n'], fsf.regstandard_nonlinear_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# High pass filter cutoff' ...
      '\nset fmri(paradigm_hp) %d\n'], fsf.paradigm_hp);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Number of lower-level copes feeding into higher-level analysis' ...
      '\nset fmri(ncopeinputs) %d\n'], fsf.ncopeinputs);
fprintf(fid,'%s', fsf_str);

nfiles = length(fsf.feat_files);
for ifile = 1:nfiles
  fsf_str = sprintf([...
	'\n# 4D AVW data or FEAT directory (%d)' ...
	'\nset feat_files(%d) "%s"\n'], ifile, ifile, fsf.feat_files{ifile});
  fprintf(fid,'%s', fsf_str);
end

if fsf.reginitial_highres_yn
  nfiles = length(fsf.initial_highres_files);
  for ifile = 1:nfiles
    fsf_str = sprintf([...
	  '\n# Session''s structural image for analysis %d' ...
	  '\nset initial_highres_files(%d) "%s"\n'], ...
	ifile, ifile, fsf.initial_highres_files{ifile});
    fprintf(fid,'%s', fsf_str);
  end
end

if fsf.reghighres_yn
  nfiles = length(fsf.highres_files);
  for ifile = 1:nfiles
    fsf_str = sprintf([...
	  '\n# Subject''s structural image for analysis %d' ...
	  '\nset highres_files(%d) "%s"\n'], ifile, ifile, fsf.highres_files{ifile});
    fprintf(fid,'%s', fsf_str);  
  end
end
  
% Here we loop over all of the variables specified in the model
nev = length(fsf.ev);
for iev = 1:nev
  ev = fsf.ev(iev);

  if fsf.version == 5.98
    fsf_str = sprintf([...
        '\n# EV %d title' ...
        '\nset fmri(evtitle%d) "%s"\n'], iev, iev, ev.name);
    fprintf(fid,'%s', fsf_str);
  end
  
  fsf_str = sprintf([...
	'\n# Basic waveform shape (EV %d)' ...
	'\n# 0 : Square' ...
	'\n# 1 : Sinusoid' ... 
	'\n# 2 : Custom (1 entry per volume)' ...
	'\n# 3 : Custom (3 column format)' ...
	'\n# 4 : Interaction' ...
	'\nset fmri(shape%d) %d\n'], iev, iev, ev.shape);
  fprintf(fid,'%s', fsf_str);
  
  fsf_str = sprintf([...
	'\n# Convolution (EV %d)' ...
	'\n# 0 : None' ...
	'\n# 1 : Gaussian' ...
	'\n# 2 : Gamma' ...
	'\n# 3 : Double-Gamma HRF' ...
	'\n# 4 : Gamma basis functions' ...
	'\n# 5 : Sine basis functions' ...
	'\n# 6 : FIR basis functions' ...
	'\nset fmri(convolve%d) %d\n'], iev, iev, ev.convolve);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
	'\n# Convolve phase (EV %d)' ...
	'\nset fmri(convolve_phase%d) %d\n'], iev, iev, ev.convolve);
  fprintf(fid,'%s', fsf_str);
  
  fsf_str = sprintf([...
	'\n# Apply temporal filtering (EV %d)' ...
	'\nset fmri(tempfilt_yn%d) %d\n'], iev, iev, ev.tempfilt_yn);
  fprintf(fid,'%s', fsf_str);
  
  fsf_str = sprintf([...
	'\n# Add temporal derivative (EV %d)' ...
	'\nset fmri(deriv_yn%d) %d\n'], iev, iev, ev.deriv_yn);
  fprintf(fid,'%s', fsf_str);
  
  switch ev.shape
    case 0
      fsf_str = sprintf([...
	    '\n# Skip (EV %d)' ...
	    '\nset fmri(skip%d) %d\n'], iev, iev, ev.skip);
      fprintf(fid,'%s', fsf_str);
      
      fsf_str = sprintf([...
	    '\n# Off (EV %d)' ...
	    '\nset fmri(off%d) %1.2f\n'], iev, iev, ev.off);
      fprintf(fid,'%s', fsf_str);
      
      fsf_str = sprintf([...
	    '\n# On (EV %d)' ...
	    '\nset fmri(on%d) %1.2f\n'], iev, iev, ev.on);
      fprintf(fid,'%s', fsf_str);
      
      fsf_str = sprintf([...
	    '\n# Phase (EV %d)' ...
	    '\nset fmri(phase%d) %d\n'], iev, iev, ev.phase);
      fprintf(fid,'%s', fsf_str);
      
      fsf_str = sprintf([...
	    '\n# Stop (EV %d)' ...
	    '\nset fmri(stop%d) %d\n'], iev, iev, ev.stop);
      fprintf(fid,'%s', fsf_str);
    case {2,3}
      fsf_str = sprintf([...
	    '\n# Custom EV file (EV %d)' ...
	    '\nset fmri(custom%d) "%s"\n'], iev, iev, ev.fname);
      fprintf(fid,'%s', fsf_str);
    otherwise
      fprintf('Not set up to handle this EV shape yet\n');
  end % switch ev.shape
  
  % Add convolution info
  switch ev.convolve
    case 1  % gaussian
      fsf_str = sprintf([...
	    '\n# Gauss sigma (EV %d)' ...
	    '\nset fmri(gausssigma%d) %1.2f\n'], iev, iev, ev.gausssigma);
      fprintf(fid,'%s', fsf_str);

      fsf_str = sprintf([...
	    '\n# Gauss delay (EV %d)' ...
	    '\nset fmri(gaussdelay%d) %1.2f\n'], iev, iev, ev.gaussdelay);
      fprintf(fid,'%s', fsf_str);
      
  end
  
  
  % Add the orthogonalization info
  if ~isempty(ev.ortho)
    ortho_vect = [0 ev.ortho];
    for jev = 0:nev
      fsf_str = sprintf([...
	    '\n# Orthogonalise EV %d wrt EV %d' ...
	    '\nset fmri(ortho%d.%d) %d'], iev, jev, iev, jev, ortho_vect(jev+1));
      fprintf(fid,'%s', fsf_str);
    end % for jev
  end
end % for iev

% Add contrast information from fsf.con structure
fsf_str = sprintf([...
      '\n# Contrast & F-tests mode' ...
      '\n# real : control real EVs' ...
      '\n# orig : control original EVs' ...
      '\nset fmri(con_mode_old) %s' ...
      '\nset fmri(con_mode) %s'], fsf.con_mode_old, fsf.con_mode);
fprintf(fid,'%s', fsf_str);

% Loop over each of the real contrasts
ncon = length(fsf.con);
nFtest = length(fsf.Ftest);

for icon = 1:ncon
  con = fsf.con(icon);
  
  mode = con.type;
  
  fsf_str = sprintf([...
	'\n# Display images for contrast_%s %d' ...
	'\nset fmri(conpic_%s.%d) %d\n'], mode, icon, mode, icon, con.conpic);
  fprintf(fid,'%s', fsf_str);

  fsf_str = sprintf([...
	'\n# Title for contrast_%s %d' ...
	'\nset fmri(conname_%s.%d) "%s"\n'], mode, icon, mode, icon, con.conname);
  fprintf(fid,'%s', fsf_str);

  for iev = 1:nev
    if length(con.con_vect) < iev, continue, end
    fsf_str = sprintf([...
	  '\n# Real contrast_%s vector %d element %d' ...
	  '\nset fmri(con_%s%d.%d) %d\n'], mode, icon, iev, mode, icon, iev, con.con_vect(iev));
    fprintf(fid,'%s', fsf_str);
  end % for iev

  for iftest = 1:length(con.ftest)
    fsf_str = sprintf([...
	  '\n# F-test %d element %d' ...
	  '\nset fmri(ftest_%s%d.%d) %d\n'], iftest, icon, mode, iftest, icon, con.ftest(iftest));
    fprintf(fid,'%s', fsf_str);
  end
end % for icon=

% Contrast masking
fsf_str = sprintf([...
      '\n# Contrast masking - use >0 instead of thresholding?' ...
      '\nset fmri(conmask_zerothresh_yn) %d\n'], fsf.conmask_zerothresh_yn);
fprintf(fid,'%s', fsf_str);

fsf_str = sprintf([...
      '\n# Do contrast masking at all?' ...
      '\nset fmri(conmask1_1) %d\n'], fsf.conmask1_1);
fprintf(fid,'%s', fsf_str);

if fsf.version == 5.98
  % Extra non-GUI options
  fsf_str = sprintf([...
      '\n#############################################################' ...
      '\n# Now options that don''t appear in the GUI\n' ...
      '\n# Alternative example_func image (not derived from input 4D dataset' ...
      '\nset fmri(alternative_example_func) %s\n' ...
      '\n# Alternative (to BETting) mask image' ...
      '\nset fmri(alternative_mask) %s\n' ...
      '\n# Initial structural space registration initialisation transform' ...
      '\nset fmri(init_initial_highres) %s\n' ...
      '\n# Structural space registration initialisation transform' ...
      '\nset fmri(init_highres) %s\n' ...
      '\n# Standard space reigstration initialisation transform' ...
      '\nset fmri(init_standard) %s\n' ...
      '\n# For full FEAT analysis: overwrite .feat output dir?' ...
      '\nset fmri(overwrite_yn) %d\n'],...
      fsf.alternative_example_func,fsf.alternative_mask,...
      fsf.init_initial_highres,fsf.init_highres,...
      fsf.init_standard,fsf.overwrite_yn);
  fprintf(fid,'%s', fsf_str);
end

fclose(fid);
return


'\n# Mask real contrast/F-test 1 with real contrast/F-test 2?'
'set fmri(conmask1_2) 0'

'\n# Mask real contrast/F-test 1 with real contrast/F-test 3?'
'set fmri(conmask1_3) 0'

'\n# Mask real contrast/F-test 1 with real contrast/F-test 4?'
'set fmri(conmask1_4) 0'

'\n# Mask real contrast/F-test 2 with real contrast/F-test 1?'
'set fmri(conmask2_1) 0'

'\n# Mask real contrast/F-test 2 with real contrast/F-test 3?'
'set fmri(conmask2_3) 0'

'\n# Mask real contrast/F-test 2 with real contrast/F-test 4?'
'set fmri(conmask2_4) 0'

'\n# Mask real contrast/F-test 3 with real contrast/F-test 1?'
'set fmri(conmask3_1) 0'

'\n# Mask real contrast/F-test 3 with real contrast/F-test 2?'
'set fmri(conmask3_2) 0'

'\n# Mask real contrast/F-test 3 with real contrast/F-test 4?'
'set fmri(conmask3_4) 0'

'\n# Mask real contrast/F-test 4 with real contrast/F-test 1?'
'set fmri(conmask4_1) 0'

'\n# Mask real contrast/F-test 4 with real contrast/F-test 2?'
'set fmri(conmask4_2) 0'

'\n# Mask real contrast/F-test 4 with real contrast/F-test 3?'
'set fmri(conmask4_3) 0'

