function makeplot = sizer2ss(data,paramstruct) 
% SIZER2SS, SIZER, version 2, for Smoothing Splines
%     Modification of SiZer2, where the underlying smoother 
%     is the cubic smoothing spline
% Inputs:
%   data        - n x 2 matrix of regression data:
%                        X's in first column,  Y's in second
%   paramstruct - a Matlab structure of input parameters
%                    Use: "help struct" and "help datatypes" to
%                         learn about these.
%                    Create one, using commands of the form:
%
%       paramstruct = struct('field1',values1,...
%                            'field2',values2,...
%                                             ) ;
%
%                          where any of the following can be used,
%                          these are optional, misspecified values
%                          revert to defaults
%
%    fields            values
%
%    iout             1  (default)  use 3 panels:  family overlay & 
%                                          surface, SiZer
%                     2  use 2 panels: family overlay, SiZer
%                     3  family overlay only
%                     4  family surface only
%                     5  SiZerSS only
%
%    imovie           1  (default)  make movie version
%                     0  make a single still plot
%                            (reset to this for iout >= 4)
%
%    icolor           1  (default)  full color version
%                     0  fully black and white version
%                     2  digital map %%added by JT 
%
%    savestr          string controlling saving of output,
%                         either a full path, or a file prefix to
%                         save in matlab's current directory
%                     unspecified:  results only appear on screen
%                     movie version (imovie = 1):
%                         add .mpg and create MPEG file
%                     static version (imovie = 0):  
%                         add .ps, and save as either
%                              color postscript (icolor = 1)
%                         or
%                              black&white postscript (when icolor = 0)
%
%    xlabelstr        String for labeling x axes (default is none)
%
%    ylabelstr        String for labeling x axes (default is none)
%
%    labelfontsize    Font Size for labels (uses Matlab default)
%                                   (15 is "fairly large")
%
%    famoltitle       Title for family overlay plot
%                           (default is 'Family Overlay, date')
%
%    famsurtitle      Title for family surface plot
%                           (default is 'SiZerSS colored Scale Space')
%
%    sizertitle       Title for SiZerSS plot
%                           (default is 'SiZerSS Map')
%    vdflag           Flag for using usual spar (0, default) or EDF (1)
%
%    titlefontsize    Font Size for titles (uses Matlab default)
%                                   (18 is "fairly large")
%
%    viewangle        Angle for Viewing of 3d surface plot, 
%                         in degrees azimuth and elevation, 
%                         recommend  [165,30] (default)
%                                or  [195,30]
%
%    ndataoverlay     0  no data plot
%                     1  (default) plot 1000 points (random choice)
%                     2  plot full data set
%                     n > 2  plot n random points
%
%    iscreenwrite     0  (default)  no screen writes
%                     1  write to screen to show progress
%
%    nrepeat          number of times to repeat movie (default = 5)
%
%    nbin             number of bins (default = 401)
%
%    minx             left end of bin range (default is min of data)
%
%    maxx             right end of bin range (default is max of data)
%
%    bpar             bin range boundary handling parameter
%                           0 - (default), move data to ends
%                           1 - truncate data outside ends
%
%    alpha            Usual "level of significance".  
%                         I.e. C.I.s have coverage probability 
%                         1 - alpha.  (0.05 when not specified)
%
%    nfspar              number of smoothing parameter's for family
%                           (default, 41 for movies,  11 for static)
%
%    fsparmin            left end of family smoothing parameter range
%                           (default = binwd / 2) 
%
%    fsparmax            right end of family smoothing parameter range
%                           (default = range)
%
%    nsspar              number of smoothing parameter's for SiZer/SiCon
%                           (default, 41 for movies,  11 for static)
%
%    ssparmin            left end of SiZerSS smoothing parameter range
%                           (default = binwd / 2) 
%
%    ssparmax            right end of SiZerSS smoothing parameter range
%                           (default = range)
%
%    ihhighlight      0 (default) don't highlight data driven h 
%                                   (in static output)
%                     1 highlight data driven h (GCV)
%
% Outputs:
%     For iout = 1,2:   graphics in current Figure
%     For iout = 3,4,5:   graphics in current axes
%     When savestr exists,
%        For imovie = 1:  MPEG file saved in 'savestr'.mpg
%        For imovie = 0:  Postscript file saved in 'savestr'.ps
%                        (color postscript for icolor = 1)
%                        (B & W postscript for icolor = 0)
%    
%
% Assumes path can find personal functions:
%  gplbinr.m
%  vec2mat.m
%  phiinv.m
%  SSfamfit.m
%  SSfamfitlim.m
%

%    Copyright (c) J. S. Marron and Jin-Ting Zhang 2001



%  First set all parameter to defaults
%
iout = 1 ;
imovie = 1 ;
icolor = 1 ;
%% change to 0 (black-white map) from 1 (color map)
edflag = 0; %% use usual spar
savestr = [] ;
ndataoverlay = 1 ;
iscreenwrite = 0 ;
nbin = min([length(data(:,1)),401]) ; %% original 401
%nbin=301;
minx = [] ;
maxx = [] ;
bpar = 0 ;
alpha = 0.05 ;
nfspar = [] ;
fsparmin = [] ;
fsparmax = [] ;
nsspar = [] ;
ssparmin = [] ;
ssparmax = [] ;
ihhighlight = 0 ;
nrepeat = 5 ;
viewangle = [165,30] ;
xlabelstr = [] ;
ylabelstr = [] ;
labelfontsize = [] ;
famoltitle = ['Family Overlay, ' date] ;
famsurtitle = ['SiZerSS colored Scale Space'] ;
sizertitle = ['SiZerSS Map'] ;
titlefontsize = [] ;





%  Now update parameters as specified,
%  by parameter structure (if it is used)
%
if nargin > 1 ;   %  then paramstruct is an argument

  if isfield(paramstruct,'iout') ;    %  then change to input value
    iout = getfield(paramstruct,'iout') ; 
  end ;

  if isfield(paramstruct,'imovie') ;    %  then change to input value
    imovie = getfield(paramstruct,'imovie') ; 
  end ;

  if isfield(paramstruct,'icolor') ;    %  then change to input value
    icolor = getfield(paramstruct,'icolor') ; 
  end ;

 if isfield(paramstruct,'edflag') ;    %  then change to input value
    edflag = getfield(paramstruct,'edflag') ; 
  end ;
  if isfield(paramstruct,'savestr') ;    %  then use input value
    savestr = getfield(paramstruct,'savestr') ; 
    if ~ischar(savestr) ;    %  then invalid input, so give warning
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from sizer2.m:     !!!') ;
      disp('!!!   Invalid savestr,           !!!') ;
      disp('!!!   using default of no save   !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      savestr = [] ;
    end ;
  end ;

  if isfield(paramstruct,'ndataoverlay') ;    %  then change to input value
    ndataoverlay = getfield(paramstruct,'ndataoverlay') ; 
  end ;

  if isfield(paramstruct,'iscreenwrite') ;    %  then change to input value
    iscreenwrite = getfield(paramstruct,'iscreenwrite') ; 
  end ;

  if isfield(paramstruct,'nbin') ;    %  then change to input value
    nbin = getfield(paramstruct,'nbin') ; 
  end ;

  if isfield(paramstruct,'minx') ;    %  then change to input value
    minx = getfield(paramstruct,'minx') ; 
  end ;

  if isfield(paramstruct,'maxx') ;    %  then change to input value
    maxx = getfield(paramstruct,'maxx') ; 
  end ;

  if isfield(paramstruct,'bpar') ;    %  then change to input value
    bpar = getfield(paramstruct,'bpar') ; 
  end ;

  if isfield(paramstruct,'alpha') ;    %  then change to input value
    alpha = getfield(paramstruct,'alpha') ; 
  end ;

  if isfield(paramstruct,'nfspar') ;    %  then change to input value
    nfspar = getfield(paramstruct,'nfspar') ; 
  end ;
  
  if isfield(paramstruct,'nsspar') ;    %  then change to input value
    nsspar = getfield(paramstruct,'nsspar') ; 
  end ;

  if isfield(paramstruct,'fsparmin') ;    %  then change to input value
    fsparmin = getfield(paramstruct,'fsparmin') ; 
  end ;

  if isfield(paramstruct,'fsparmax') ;    %  then change to input value
    fsparmax = getfield(paramstruct,'fsparmax') ; 
  end ;

  

  if isfield(paramstruct,'ssparmin') ;    %  then change to input value
    ssparmin = getfield(paramstruct,'ssparmin') ; 
  end ;

  if isfield(paramstruct,'ssparmax') ;    %  then change to input value
    ssparmax = getfield(paramstruct,'ssparmax') ; 
  end ;

  if isfield(paramstruct,'ihhighlight') ;    %  then change to input value
    ihhighlight = getfield(paramstruct,'ihhighlight') ; 
  end ;

  if isfield(paramstruct,'viewangle') ;    %  then change to input value
    viewangle = getfield(paramstruct,'viewangle') ; 
  end ;

  if isfield(paramstruct,'nrepeat') ;    %  then change to input value
    nrepeat = getfield(paramstruct,'nrepeat') ; 
  end ;

  if isfield(paramstruct,'xlabelstr') ;    %  then change to input value
    xlabelstr = getfield(paramstruct,'xlabelstr') ; 
  end ;

  if isfield(paramstruct,'ylabelstr') ;    %  then change to input value
    ylabelstr = getfield(paramstruct,'ylabelstr') ; 
  end ;

  if isfield(paramstruct,'labelfontsize') ;    %  then change to input value
    labelfontsize = getfield(paramstruct,'labelfontsize') ; 
  end ;

  if isfield(paramstruct,'famoltitle') ;    %  then change to input value
    famoltitle = getfield(paramstruct,'famoltitle') ; 
  end ;

  if isfield(paramstruct,'famsurtitle') ;    %  then change to input value
    famsurtitle = getfield(paramstruct,'famsurtitle') ; 
  end ;

  if isfield(paramstruct,'sizertitle') ;    %  then change to input value
    sizertitle = getfield(paramstruct,'sizertitle') ; 
  end ;

  if isfield(paramstruct,'titlefontsize') ;    %  then change to input value
    titlefontsize = getfield(paramstruct,'titlefontsize') ; 
  end ;


end ;  %  of resetting of input parameters




%  Get Matlab version
%
verstr = version ;
vernum = str2num(verstr(1:3)) ;



%  Turn iout into control parameters
%
if iout == 1 ;
  viplot = [1; 1; 1] ;
         %  indicators for family overlay, family surface, SiZer, SiCon
elseif iout == 2 ;
  viplot = [1; 0; 1] ;
elseif iout == 3 ;
  viplot = [1; 0; 0] ;
elseif iout == 4 ;
  viplot = [0; 1; 0] ;
elseif iout == 5 ;
  viplot = [0; 0; 1] ;
end ;
nplot = sum(viplot) ;

if  nplot == 1  &  imovie == 1  ;
          %  then reset imovie, since it focuses on figures, not axes 
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Warning from sizer2ss.m:   !!!') ;
  disp('!!!   resetting imovie to 0      !!!') ;
  disp('!!!   (use iout < 4 for movie)   !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  imovie = 0 ;
end ;



%  set movie parameters
%
nrif = 4 ;
          %  number of times to Repeat Individual Movie Frames
          %    (controls speed of .mpeg version)
movopt = 1 ;
%  Movie options, currently using same as Cornean
%
%  This is MATLAB "help" to use for choosing parameters:
%	MPGWRITE(M, map, 'filename', options) Encodes M in MPEG format
% 	using the specified colormap and writes the result to the
% 	specified file.  The options argument is an optional vector of
% 	8 or fewer options where each value has the following meaning:
% 	1. REPEAT:
% 		An integer number of times to repeat the movie
% 		(default is 1).
% 	2. P-SEARCH ALGORITHM:
% 		0 = logarithmic	(fastest, default value)
% 		1 = subsample
% 		2 = exhaustive	(better, but slow)
% 	3. B-SEARCH ALGORITHM:
% 		0 = simple	(fastest)
% 		1 = cross2	(slightly slower, default value)
% 		2 = exhaustive	(very slow)
% 	4. REFERENCE FRAME:
% 		0 = original	(faster, default)
% 		1 = decoded	(slower, but results in better quality)
% 	5. RANGE IN PIXELS:
% 		An integer search radius.  Default is 10.
% 	6. I-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 8.
% 	7. P-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 10.
% 	8. B-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 25.




xdat = data(:,1) ;
ydat = data(:,2) ;
mdat = [xdat,ydat] ;



%  check not all data points are same
%
if std(xdat) == 0 ;     %  then all data points are same, warn and quit
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from sizer2ss.m:   !!!') ;
  disp('!!!   All x values are same,   !!!') ;
  disp('!!!   Terminating Execution    !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  return ;
end ;




%  Determine data ranges
%
if isempty(minx) ;
  minx = min(xdat) ;
end ;

if isempty(maxx) ;
  maxx = max(xdat) ;
end ;

ndat = length(xdat) ;




%  Determine smoothing parameter ranges
%
if isempty(nfspar) ;
  if imovie == 0 ;
    nfspar = 11 ; %% original 11
  else ;
    nfspar = 41 ;
  end ;
end ;



if isempty(fsparmin) ;
  if vernum >= 6.0 ;
    %fsparmin = 10^fzero('SSfamfitlim',1,[],nbin,0.8*nbin,[minx,maxx]); 
  else ;
    %fsparmin = 10^fzero('SSfamfitlim',1,[],[],nbin,0.8*nbin,[minx,maxx]); 
  end ;
  
end ;


if isempty(fsparmax) ;
  if vernum >= 6.0 ;
    %fsparmax = 10^fzero('SSfamfitlim',1,[],nbin,2.2,[minx,maxx]); 
  else ;
    %fsparmax = 10^fzero('SSfamfitlim',1,[],[],nbin,2.2,[minx,maxx]); 
  end ;
  
end ;
temp=sparband(mdat(:,1),nfspar,nbin);
fsparmin=temp(1); fsparmax=temp(nfspar);




if nfspar > 1 ;
  if fsparmin >= fsparmax ;
    disp('!!!  Warning from sizer2ss:                 !!!') ;
    disp('!!!  nfspar > 1, and fsparmin >= fsparmax      !!!') ;
    disp('!!!  will reset nfspar to 1, & use fsparmax   !!!') ;
    nfspar = 1 ;
    fvspar = fsparmax ;
  else ;
    fvspar = logspace(log10(fsparmin),log10(fsparmax),nfspar) ;
  end ;
else ;
  if fsparmin == fsparmax ;
    fvspar = fsparmax ;
  else ;
    disp('!!!  Warning from sizer2ss:                     !!!') ;
    disp('!!!  nfspar = 1, and fsparmin, fsparmax different  !!!') ;
    disp('!!!  will use fsparmax as single h               !!!') ;
    fvspar = fsparmax ;
  end ;
end ;


if isempty(nsspar) ;
  if imovie == 0 ;
    nsspar = 11 ;  %% original 11
  else ;
    nsspar = 41 ;%% original 41
  end ;
end ;

if isempty(ssparmin) ;
  if vernum >= 6.0 ;
    %ssparmin = 10^fzero('SSfamfitlim',1,[],nbin,0.8*nbin,[minx,maxx]) ;
  else ;
    %ssparmin = 10^fzero('SSfamfitlim',1,[],[],nbin,0.8*nbin,[minx,maxx]) ;
  end ;
end ;


if isempty(ssparmax) ;
  if vernum >= 6.0 ;
    %ssparmax = 10^fzero('SSfamfitlim',1,[],nbin,2.2,[minx,maxx]) ;
  else ;
    %ssparmax = 10^fzero('SSfamfitlim',1,[],[],nbin,2.2,[minx,maxx]) ;
  end ;
end ;

temp=sparband(mdat(:,1),nsspar,nbin);
ssparmin=temp(1); ssparmax=temp(nsspar);

if nsspar > 1 ;
  if ssparmin >= ssparmax ;
    disp('!!!  Warning from sizer2ss:                 !!!') ;
    disp('!!!  nsspar > 1, and ssparmin >= ssparmax      !!!') ;
    disp('!!!  will reset nsspar to 1, & use ssparmax   !!!') ;
    nsspar = 1 ;
    svh = ssparmax ;
  else ;
    svh = logspace(log10(ssparmin),log10(ssparmax),nsspar) ;
  end ;
else ;
  if ssparmin == ssparmax ;
    svh = ssparmax ;
  else ;
    disp('!!!  Warning from sizer2ss:                     !!!') ;
    disp('!!!  nsspar = 1, and ssparmin, ssparmax different  !!!') ;
    disp('!!!  will use ssparmax as single h               !!!') ;
    svh = ssparmax ;
  end ;
end ;



%  Set up colormap
%
if icolor ~= 0 ;     %  Then go for nice colors in sizer and sicon

  %  Set up colorful color map
  cocomap = [0,    0,   1; ...
            .35, .35, .35; ...
            .5,    0,  .5; ...
             1,    0,   0; ...
             1,   .5,   0; ...
           .35,  .35, .35; ...
             0,    1,   0; ...
             0,    1,   1] ; 
  colormap(cocomap) ;

  famcolorstr = 'b' ;
          %  Plot most curves in blue
  highlightcolorstr = 'w' ;
          %  Plot highlighted curve in white
  datcolorstr = 'g' ;
          %  Plot overlaid data in green

else ;     %  Then use gray scale maps everywhere

  %  Set up gray level color map
  comap = [.2, .2, .2; ...
           .35, .35, .35; ...
           .5, .5, .5; ...
           .8, .8, .8] ;
  colormap(comap) ;

  famcolorstr = 'k' ;
%%orig 'w'
          %  Plot most curves in white
  highlightcolorstr = 'k' ;
%% orig 'w'
          %  Plot highlighted curve in white
  datcolorstr = 'k' ;
%% orig 'w'
          %  Plot overlaid data in white

end ;





%  Calculate family of smooths (if needed)
%
if  viplot(1) == 1  |  viplot(2) == 1  ;    %  Then will show a family plot

  if iscreenwrite == 1 ;
    disp('    SiZerSS  Working on family') ;
  end ;



  [fxgrid, mfam] = SSfamfit(mdat,fvspar,nbin,[minx,maxx]) ;
  fbottom = min(min(mfam)) ;

  ftop = max(max(mfam)) ;
  frange = ftop - fbottom ;
  fbottom = fbottom - 0.05 * frange ;
  ftop = ftop + 0.05 * frange ;


end ;    %  of family plot construction if block






%  Calculate data driven bandwidths for highlighting (if needed)
%
if  imovie == 0   & ...
    ihhighlight == 1  ;

  [xfit,fhat,fsig,dhat,dsig,ledf,spargcvcurve] = ...
                          SSfamfit(mdat,fvspar,nbin,[minx,maxx]) ;

  spargcv = 10^gpminr(log10(fvspar)',log10(spargcvcurve(:,2))) ;

  if spargcv == fsparmin ;    %  if data based h below range
    disp('!!!   Warning from sizer2ss: databased h below this range   !!!') ;
    spargcvflag = 0 ;
  elseif spargcv == fsparmax ;    %  if databased h above this range
    disp('!!!   Warning from sizer2ss: databased h above this range   !!!') ;
    spargcvflag = 0 ;
  else ;     %  if data based in range, then get its nearest index
    [temp, hcind] = min(abs(log10(spargcv) - log10(fvspar))) ;
    spargcvflag = 1 ; 
  end ;

else ;

  spargcvflag = 0 ;

end ;




%  Calculate SiZerSS (if needed)
%
if  viplot(2) == 1  |  viplot(3) == 1  ; 

  if iscreenwrite == 1 ;
    disp('    SiZerSS  Working on SiZer') ;
  end ;

  disp(['# of spar for SiZerSS=', num2str(nsspar)])
  
  if edflag==0, %% use usual spar
    sizermap= gpsz1ss(mdat,[minx,maxx,nbin],[ssparmin,ssparmax,nsspar],alpha) ;
  elseif edflag==1, %% use EDF
    [sizermap,temp,vedf] = gpsz1ss(mdat,[minx,maxx,nbin],[ssparmin,ssparmax,nsspar],alpha) ;
  end
end ;



%  Do family overlay graphics (if needed) 
%
if viplot(1) == 1 ;

  if nplot > 1 ;    %  then doing multiple graphics
     
    fighand = gcf ;
    if nplot == 3 ;
      famolh = subplot(2,2,1) ;
    else ;
      famolh = subplot(nplot,1,1) ;
    end ;
  end ;


  %  Plot family
  %
  vfamh = plot(fxgrid,mfam,famcolorstr) ;
          %  vfamh = vector of family curve handles     
    vax = [minx,maxx,fbottom,ftop] ;
    axis(vax) ;
    th = title(famoltitle) ;
    if ~isempty(titlefontsize) ;
      set(th,'FontSize',titlefontsize) ;
    end ;
    lxh = xlabel(xlabelstr) ;
    lyh = ylabel(ylabelstr) ;
    if ~isempty(labelfontsize) ;
      set(lxh,'FontSize',labelfontsize) ;
      set(lyh,'FontSize',labelfontsize) ;
    end ;


  %  Highlight data driven curve (if needed)
  %
  if spargcvflag == 1 ;
    set(vfamh(hcind),'LineWidth',2) ;
          %  use fatter line width for data based choice
    set(vfamh(hcind),'Color',highlightcolorstr) ;
          %  use red color for one in middle
  end ;


  %  Overlay data (if needed)
  %
  if ndataoverlay > 0 ;    %  then overlay a plot of the raw data

    seed = 34593874 ;
    rand('seed',seed) ;

    if ndataoverlay == 1 ;
      ndo = min(ndat,1000) ;
    elseif ndataoverlay == 2 ;
      ndo = ndat ;
    else ;
      ndo = min(ndat,ndataoverlay) ;
    end ;


    if ndo < ndat ;    %  then need to subsample
      [temp,randperm] = sort(rand(ndat,1)) ;
            %  randperm is a random permutation of 1,2,...,ndat
      vkeep = randperm(1:ndo) ;
            %  indices of elements to keep for display
      mdatol = mdat(vkeep,:) ;
    else ;    %  overlay full data set
      mdatol = mdat ;
    end ;


    %  overlay selected data
    %
    hold on ;
      plot(mdatol(:,1),mdatol(:,2),[datcolorstr '.']) ;
    hold off ;
  end ;
  %%colorbar
end ;








%  Do family surface graphics (if needed) 
%
if viplot(2) == 1 ;

  if nsspar == nfspar ;    % then go ahead with surface graphics

    if nplot > 1 ;    %  then doing multiple graphics
      famsurh = subplot(2,2,2) ;
            %  this only appears in the 3 panel plot
    end ;

    %  First do decimation if needed
    %
    if nbin > 41 ;

      decimfact = ceil((nbin - 1) / 40) ;
            %  factor to decimate by
      keepflag = (1:decimfact:nbin)' ;

      mfamdecim = mfam(keepflag,:) ;
      fxgriddecim = fxgrid(keepflag) ;
      sizermapdecim = sizermap(:,keepflag) ;

    else ;

      mfamdecim = mfam ;
      fxgriddecim = fxgrid ;
      sizermapdecim = sizermap ;

    end ;


    l10hdecim = log10(fvspar) ;
    if nfspar > 21 ;

      decimfact = ceil((nfspar - 1) / 20) ;
            %  factor to decimate by
      keepflag = (1:decimfact:nfspar)' ;

      mfamdecim = mfamdecim(:,keepflag) ;
      sizermapdecim = sizermapdecim(keepflag,:) ;
      l10hdecim = l10hdecim(keepflag) ;

    end ;


    nrowsmd = size(sizermapdecim,1) ;
    ncolsmd = size(sizermapdecim,2) ;
    sizermapdecim = sizermapdecim(2:nrowsmd,2:ncolsmd) ;
          %  cut off one row and one column, to give number
          %  of panels

    %  Make surface plot
    %
    surf(fxgriddecim,l10hdecim,mfamdecim',sizermapdecim) ;
          %  surface plot, using SiZerSS colors
      th = title(famsurtitle) ;
      if ~isempty(titlefontsize) ;
        set(th,'FontSize',titlefontsize) ;
      end ;
      lxh = xlabel(xlabelstr) ;
      lyh = ylabel('log10(\lambda)') ;
      lzh = zlabel(ylabelstr) ;
      if ~isempty(labelfontsize) ;
        set(lxh,'FontSize',labelfontsize) ;
        set(lyh,'FontSize',labelfontsize) ;
        set(lzh,'FontSize',labelfontsize) ;
      end ;

      vax = [minx,maxx,log10(ssparmin),log10(ssparmax),fbottom,ftop] ;
      axis(vax) ;
      set(gca,'Xgrid','off') ;
      set(gca,'Ygrid','off') ;
      set(gca,'Zgrid','off') ;
      set(gca,'Xdir','reverse')
      view(viewangle) ; 

      caxis([1, 8]) ;
      if icolor ~= 0 ;     %  Then go for nice colors in sizer and sicom
        colormap(cocomap) ;  
      else ;     %  Then use gray scale maps everywhere
        colormap(comap) ;  
      end ;
    else ;

      disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;
      disp(['!!!   Error from sizer2ss.m:    !!!']) ;
      disp(['!!!   For surface plot,         !!!']) ;
      disp(['!!!   must have nsspar = nfspar   !!!']) ;
      disp(['!!!   Terminating execution     !!!']) ;
      disp(['!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!']) ;
      return ;

    end ;


  end ;    %  of family surface plot construction if block






%  Plot SiZerSS (if needed)
%
if viplot(3) == 1 ;    %  Then will show a SiZerSS map

  if nplot > 1 ;    %  then doing multiple graphics
    if nplot == 3 ;
      sizerh = subplot(2,2,3) ;
    else ;
      sizerh = subplot(nplot,1,2) ;
    end ;
  end ;

    
   if edflag==0, %% use usual spar
       image([minx,maxx],[log10(ssparmin),log10(ssparmax)],sizermap) ;
   elseif edflag==1, %% use effective degrees of freedom   
       image([minx,maxx],vedf,sizermap) ;
 
   end
        
    set(gca,'YDir','normal') ;
    th = title(sizertitle) ;
    if ~isempty(titlefontsize) ;
      set(th,'FontSize',titlefontsize) ;
    end ;
    lxh = xlabel(''); %%xlabel(xlabelstr) ;%%removed by zjt nov 18,2004
    
     if edflag==0,
        lyh = ylabel('log10(\lambda)') ;
     elseif edflag==1,   
        lyh = ylabel('EDF');
     end   
     
    if ~isempty(labelfontsize) ;
      set(lxh,'FontSize',labelfontsize) ;
      set(lyh,'FontSize',labelfontsize) ;
    end ;
          %  next add marking lines
    if icolor ~= 0 ;     %  Then go for nice colors in sizer and sicom
      colormap(cocomap) ;  
    else ;     %  Then use gray scale maps everywhere
      colormap(comap) ;  
    end ;
    hold on ;
      if spargcvflag == 1 ;
        plot([minx; maxx], ones(2,1)*log10(spargcv),['-' highlightcolorstr]) ;
      end ;
    hold off ;
    
if icolor==0,    
   colorbar('horiz');%%added by ZJT, 13/11/2001
end
   
end ;    %  of SiZerSS plot construction if block

%  Do main movie construction
%
if imovie == 1 ;

  if iscreenwrite == 1 ;
    disp('    SiZerSS  Working on Movie') ;
  end ;



  %  Do highlight on family plot (if needed)
  %
  if viplot(1) == 1 ;
    axes(famolh) ;
    hold on ;
      famollineh = plot(fxgrid,mfam(:,1),highlightcolorstr) ;
      set(famollineh,'LineWidth',2) ;
          %  use fatter line width for data based choice
    hold off ;
  end ;

  %  Do highlight on family surface (if needed)
  %
  if viplot(2) == 1 ;
    axes(famsurh) ;
    hold on ;
      famsurlineh = plot3(fxgrid,log10(fvspar(1))*ones(nbin,1), ...
                                   mfam(:,1),highlightcolorstr) ;
      set(famsurlineh,'LineWidth',2) ;
          %  use fatter line width for data based choice
    hold off ;
  end ;

  %  Add highlight to SiZerSS map (if needed)
  %
  if viplot(3) == 1 ;
    axes(sizerh) ;
          %  make SiZerSS axes current
    hold on ;
      sizerlineh = plot([minx; maxx], ones(2,1)*log10(fvspar(1)),['-',highlightcolorstr]) ;
    hold off ;
  end ;



  mmovie(:,1) = getframe(fighand) ;


  %  Loop through and change
  for iframe = 2:nfspar ;

    h = fvspar(iframe) ;

    if viplot(1) == 1 ;
      set(famollineh,'YData',mfam(:,iframe)) ;
    end ;


    if viplot(2) == 1 ;
      set(famsurlineh,'YData',log10(fvspar(iframe))*ones(nbin,1)) ;
      set(famsurlineh,'ZData',mfam(:,iframe)) ;
    end ;


    if viplot(3) == 1 ;
      set(sizerlineh,'YData',ones(2,1)*log10(fvspar(1))) ;
% org. \lambda
    end ;


    mmovie(:,iframe) = getframe(fighand) ;

  end ;    %  of iframe loop, to make movie



  %  Reorder frames, and replay
  %
  vmorder = [(1:nfspar),((nfspar-1):-1:2)] ;
  vmorder = vec2mat(vmorder,nrif) ;
  vmorder = reshape(vmorder,1,size(vmorder,2)*nrif) ;
          %  replicate frame number, NRIF times

  mmovie = mmovie(:,vmorder') ;
          %  reorder frames, to loop back to beginning

  movie(fighand,mmovie,nrepeat) ;
          %  Play movie on screen

end ; %% end for imovie==1





%  Save results (if needed)
%
if ~isempty(savestr) ;     %  then save results

  if iscreenwrite == 1 ;
    disp('    SiZerSS  saving results') ;
  end ;


  if imovie == 0 ;     %  then save as postscript file

    if  nplot == 1  |  nplot == 3 ;
      orient landscape ;
    else ;
      orient tall ;
    end ;

    if icolor ~= 0 ;     %  then make color postscript
      eval(['print -dpsc ' savestr '.ps']) ;
    else ;                %  then make black and white
      eval(['print -dps ' savestr '.ps']) ;
    end ;

  elseif imovie == 1 ;    %  then save as mpeg file

    if icolor ~= 0 ;     %  then save in full color
      mpgwrite(mmovie,cocomap,savestr,movopt) ;
    else ;                %  then save in black and white
      mpgwrite(mmovie,comap,savestr,movopt) ;
    end ;

  else ;

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from sizer2ss.m:               !!!') ;
    disp('!!!   Invalid value of imovie,             !!!') ;
    disp('!!!   Terminating Execution without save   !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    return ;

  end ;


  if iscreenwrite == 1 ;
    disp('    SiZerSS  finished save') ;
    disp('  ') ;
  end ;

end ;