disp('Running MATLAB script file ScaSpaEg2.m') ;
%
%    Scale based Stationarity examples 2
%    Simulaiton study to check size for various simul types


%addpath \matlab\stevesmo ;
%addpath \matlab\stevefun ;
addpath \matlabR12\work\SteveFun -end ;
addpath \matlabR12\work\SteveSmo -end ;


ipart = 2 ;    %  1 - Make several images (for different gaps)
               %  2 - Make movie through gaps
               

%visimul = [0 1 2 3] ;
visimul = [4] ;
                %  0 - no multiple comparison adjustment
                %  1 - adjust for multiple comparisons
                %  2 - Bonferroni adjustment
                %  3 - mixture
                %  4 - global simultaneous
               
%vidat = [1 2 3] ;
vidat = [1] ;
              %  1 - no signal
              %  2 - sine wave
              %  3 - doppler
              

%vinoise = [1 2 3 4 5 6 7] ;
vinoise = [10] ;
                %  1 - no noise
                %  2 - lo Gaussian Noise, sig = 0.1
                %  3 - hi Gaussian Noise, sig = 1 
                %  4 - AR(1), 0.99, var = 0.1
                %  5 - AR(1), 0.99, var = 1
                %  6 - sin varying AR(1)
                %  7 - AR(1), 0.75, var = 0.1
                %  8 - N(0,sin varying), var = 0.1
                %  9 - sin varying mixture of heavy tail and bimodal
                %  10 - Cauchy, quartiles at +-0.1


vitest = [2 1] ;
               %  1 - means and t-test
               %  2 - Kolmogorov-Smirnov
               %  3 - variances and F-test
               %  4 - mean removed Kolmogorov-Smirnov
               %  5 - standardized Kolmogorov-Smirnov
               %  6 - median removed Kolmogorov-Smirnov
               %  7 - median and MAD st'dized Kolmogorov-Smirnov
               %  8 - Quartile st'dized Kolmogorov-Smirnov



%  Set common parameters
%  
n = 10000 ;
left = 1 ;
right = n ;
%vgap = [1; 5; 10; 30; 100; 300] ;
%vgap = [1; 15; 200; 3000] ;
%vgap = [1 4000] ;
%vgap = [1] ;
vgap = logspace(log10(1),log10(3000),10) ;
ngap = length(vgap) ;
mingap = min(vgap) ;
nwid = 20 ;
%nwid = 40 ;
          %  number of widths to show at each scale
%nloc = 50 ;
%nloc = 400 ;
nloc = 100 ;
          %  number of locations to consider, for each scale
alpha = 0.05 ;
nminew = 5 ;
          %  minimum number in each window
nrepeat = 3 ;
          %  number of times to repeat movie (while showing)
fps = 2 ;
          %  frames per second, for movie (both show and save)
          
%  Set Colormap
%
comap = [.35, .35, .35; ...
         0,    0,   1; ...
         1,    0,   0] ; 
colormap(comap) ;

figh = figure(1) ;
clf ;



%  loop through simul levels
%
for isimul = visimul ;

  randn('seed',44616168) ;
  %randn('seed',34029874) ;
  %randn('seed',91237236) ;
  %randn('seed',91237345) ;
  %randn('seed',91287345) ;
  rand('seed',29837439) ;
      %  only used for inoise = 9 
  if isimul == 0 ;   
          simultit = 'not simultaneous' ;      
	elseif isimul == 1 ; 
          simultit = 'simultaneous' ;      
	elseif isimul == 2 ; 
          simultit = 'Bonferroni' ;    
  elseif isimul == 3 ; 
          simultit = 'mixture' ;  
  elseif isimul == 4 ; 
          simultit = 'global' ;  
          %mglobal = 2 * nloc - 1 ;
          mglobal = 2 * min((right - left)/nminew,nloc) - 1 ;
          salphaglobal = 1 - (1 - alpha)^(1/mglobal) ; 
  end;
  
  

	%  loop through tests
	%
	for itest = vitest ;
	
	
	
		if itest == 1 ;    %  1 - means and t-test
          testtit = 'Mean t-test' ;      
		elseif itest == 2 ;    %  2 - Kolmogorov-Smirnov
          testtit = 'Kolm-Smir' ;      
		elseif itest == 3 ;    %  3 - variances and F-test
            testtit = 'Var. F-test' ;      
		elseif itest == 4 ;    %  4 - mean removed Kolmogorov-Smirnov
            testtit = 'Mean R. K-S' ;      
		elseif itest == 5 ;    %  5 - standardized Kolmogorov-Smirnov
            testtit = 'Std''zed K-S' ;  
    elseif itest == 6;    %  6 - median removed Kolmogorov-Smirnov
            testtit = 'Median R. K-S' ;    
		elseif itest == 7 ;    %  7 - median and MAD st'dized Kolmogorov-Smirnov
            testtit = 'MAD S''d K-S' ;  
		elseif itest == 8 ;    %  8 - Quartile st'dized Kolmogorov-Smirnov
            testtit = 'Q S''d K-S' ;  
		end ;    %  end of itest if-block
		
	
		
		
		%  loop through data sets
		%
		for idat = vidat ;
		
		
		
			%  Generate Data
			%
			tgrid = linspace(left,right,n)' ;
			
			
			%  Generate underlying signal
			%
			if idat == 1 ;
			
              dattit = 'No Signal' ;
			
              vsignal = zeros(n,1) ;
			
			elseif idat == 2 ;
			
              dattit = 'Sine Wave' ;
			
              vsignal = sin(2 * pi * tgrid / 2000) ;
			
			elseif idat == 3 ;
			
              dattit = 'Doppler' ;
			
              vsignal = sin(2 * pi * (tgrid / n).^3 * n / 300) ;
			
			end ;    %  of idat if-block
			
			
			
			%  loop through noise levels
			%
			for inoise = vinoise ;
			
			
				
				%  Add noise
				%
				if inoise == 1 ;    %  don't add noise
				
          noisetit = 'no noise' ;
          
          vx = vsignal;
                  
				elseif inoise == 2 ;    %  lo Gaussian Noise, sig = 0.1
				
          noisetit = 'N, s=0.1 noise' ;
          
          sig = 0.1 ;
          vx = vsignal + sig * randn(n,1) ;

          
				elseif inoise == 3 ;    %  hi Gaussian Noise, sig = 1
				
          noisetit = 'N, s=1 noise' ;
          
          sig = 1 ;
          vx = vsignal + sig * randn(n,1) ;

          
				elseif inoise == 4 ;    %  AR(1), 0.99, var = 0.1
				
          noisetit = '0.99 AR(1), v0.1' ;
          
          phi = 0.99 ;
          sig0 = 0.1 ;
              %  marginal standard deviation
          sig = sig0 * sqrt(1 - phi^2) ;
              %  standard deviation of noise terms

          vnoise = sig * randn(n-1,1) ;
          
          varn = zeros(n,1) ;
          varn(1) = sig0 * randn(1,1) ;
          
          for i=2:n ;
            varn(i) = phi * varn(i-1) + vnoise(i-1) ;
          end ;

          vx = vsignal + varn ;

          
				elseif inoise == 5 ;    %  AR(1), 0.99, var = 1
				
          noisetit = '0.99 AR(1), v1' ;
          
          phi = 0.99 ;
          sig0 = 1 ;
              %  marginal standard deviation
          sig = sig0 * sqrt(1 - phi^2) ;
              %  standard deviation of noise terms

          vnoise = sig * randn(n-1,1) ;
          
          varn = zeros(n,1) ;
          varn(1) = sig0 * randn(1,1) ;
          
          for i=2:n ;
            varn(i) = phi * varn(i-1) + vnoise(i-1) ;
          end ;

          vx = vsignal + varn ;

          
				elseif inoise == 6 ;    %  sin varying AR(1)
				
          noisetit = 'sin varying AR(1)' ;
          
          %vphi = (1 - cos(2 * pi * tgrid / 5000)) / 2 ;
          vphi = sqrt((1 - cos(2 * pi * tgrid / 5000)) / 2) ;
          sig0 = 0.1 ;
              %  marginal standard deviation

          vnoise = sig0 * randn(n-1,1) ;
          
          varn = zeros(n,1) ;
          varn(1) = sig0 * randn(1,1) ;
          
          for i=2:n ;
            varn(i) = vphi(i) * varn(i-1) + vnoise(i-1) ;
          end ;

          varn = 0.1 * varn / std(varn) ;
              %  put on same time scale as "lo noise"

          vx = vsignal + varn ;
          
                    
				elseif inoise == 7 ;    %  AR(1), 0.75, var = 0.1
				
          noisetit = '0.75 AR(1), v0.1' ;
          
          phi = 0.75 ;
          sig0 = 0.1 ;
              %  marginal standard deviation
          sig = sig0 * sqrt(1 - phi^2) ;
              %  standard deviation of noise terms

          vnoise = sig * randn(n-1,1) ;
          
          varn = zeros(n,1) ;
          varn(1) = sig0 * randn(1,1) ;
          
          for i=2:n ;
            varn(i) = phi * varn(i-1) + vnoise(i-1) ;
          end ;

          vx = vsignal + varn ;
          
          
				elseif inoise == 8 ;    %  N(0,sin varying), var = 0.1

          noisetit = 'N, sin var.' ;
          
          vsig = 0.2 * (1 - cos(2 * pi * tgrid / 500)) / 2 ;
          vx = vsignal + vsig .* randn(n,1) ;
				
				
				elseif inoise == 9 ;    %  sin varying mixture of heavy tail and bimodal
				
          noisetit = 'sin mixed dist' ;

          vnoise = randn(n,1) ;
              %  N(0,1)
          vunif = rand(n,2) ;
              %  Uniform(0,1)
          vprobclass = (1 - cos(2 * pi * tgrid / 500)) / 2 ;
              %  Prob of choosing class 2
				  classflag = vunif(:,1) > vprobclass ;
				      %  1 for class 1
				      %  0 for class 2
				  mixflag = vunif(:,2) > 0.5 ;
				      %  Uniform on {0,1}
				  
				  sig1 = 0.2 ;
				  mu1 = sqrt(1 - sig1^2) ;
				      %  ensures mean = 0, and var = 1
				  
				  sig21 = 0.2 ;
				      %  smaller variance similar to above
				  sig22 = sqrt(2 - sig21^2) ;
				      %  ensures mean = 0, and var = 1
				  
				  vmu = classflag .* (mu1 * (2 * mixflag - 1)) ; 
				      %  0's for class 2
				      %  +- mu1 for class 1

				  vsig = classflag .* (sig1) + ...
				            (1 - classflag) .* (mixflag .* sig21 + (1 - mixflag) .* sig22) ;
				      %  sig1 for class 1
				      %  mixture of sig21 and sig22 for class 2
				      
				  vx = vsignal + vmu + vsig .* vnoise ;

				  
 				elseif inoise == 10 ;    %  Cauchy, quartiles at +-0.1
				
          noisetit = 'Cauchy, Q13+-0.1' ;
          
          q3 = 0.1 ;
          
          vx = vsignal + q3 * tan(pi * randn(n,1)) ;

				  
        end;    %  of noise addition if-block



				
				
        %  loop through gaps
        %
        clear moviestruct ;
        for igap = 1:ngap ;

          gap = vgap(igap) ;
          go2 = gap / 2 ;
                %  gap over 2


          %  loop through wid's
          %
          mout = zeros(nwid,nloc) ;

%for the movie need to adjust these lines (the locations 
%and widths are different across gaps)
          if ipart == 1 ;    %  just making static plot, so use full size
          
            minwid = log10(nminew) ;
            maxwid = log10((right - left - gap) / 2) ;
            vwid = logspace(minwid,maxwid,nwid+1)' ;
            vwid = vwid(1:nwid) ;
                  %  cut off last one, since doesn't add to plot
            maxwid = log10(vwid(nwid)) ;
                  %  readjust maxwid
                
            minloc = (left + nminew + go2) ;
            maxloc = (right - nminew - go2) ;

          elseif ipart == 2 ;    %  just making static plot, so use full size

            minwid = log10(nminew) ;
            maxwid = log10((right - left - mingap) / 2) ;
            vwid = logspace(minwid,maxwid,nwid+1)' ;
            vwid = vwid(1:nwid) ;
                  %  cut off last one, since doesn't add to plot
            maxwid = log10(vwid(nwid)) ;
                  %  readjust maxwid
                
            minloc = (left + nminew + mingap / 2) ;
            maxloc = (right - nminew - mingap / 2) ;
          
          end ;
            
            
          vloc = linspace(minloc,maxloc,nloc) ;
          for iwid = 1:nwid ;

            wid = vwid(iwid) ;

            if iwid == 1 ;
              disp(['  Working on wid ' num2str(iwid) ' = ' ...
                             num2str(wid) ...
                             '    for ' simultit ', ' testtit ', ' ...
                             dattit ', ' noisetit ', gap ' num2str(gap)]) ;
            end ;
                             
            locflag = left + wid + go2 <= vloc ;
            locflag = locflag & (vloc <= right - wid - go2) ;
            
            vindloc = (1:nloc) ;
            vindloc = vindloc(locflag) ;
                
	
            if isimul == 0 ;    %  don't bother about simultaneous inference
              salpha = alpha ;
            elseif isimul == 1 ;
              %m = (right - left) / (2 * wid) ;
              m = (right - left) / wid ;
              m = min(length(vindloc),m+1) ;
              salpha = 1 - (1 - alpha)^(1/m) ;

            elseif isimul == 2 ;
              salpha = alpha / length(vindloc) ;
     
            elseif isimul == 3 ;
              m = (right - left) / (2*wid) ;
              m = min(length(vindloc),m) ;
              salpha = 1 - (1 - (alpha/nwid))^(1/m) ;

            elseif isimul == 4 ;
              salpha = salphaglobal ;

            end ;
	
                
                
            for iloc = vindloc ;

              loc = vloc(iloc) ;
                
              flagleft = (loc - go2 - wid <= tgrid) ;
              flagleft = flagleft & (tgrid <= loc - go2) ;
                  %  ones where data inside left group

              flagright = (loc + go2 <= tgrid) ;
              flagright = flagright & (tgrid <= loc + go2 + wid) ;
                  %  ones where data inside right group

              dataleft = vx(flagleft) ;
              dataright = vx(flagright) ;


              if itest == 1 ;    %  means and t-test

                meanleft = mean(dataleft) ;
                meanright = mean(dataright) ;
                varleft = var(dataleft) ;
                varright = var(dataright) ;
                
                if meanleft == meanright ;                      

                  h = 0 ;
                      %  no significant difference

                else ;
                  
                  if varright == 0 ;
                    if varleft == 0 ;
                      h = 1 ;
                         %  significant difference

                    else ;    %  do 1 sample t-test
                      h = ttest(dataleft,meanright,salpha) ;
                         %  1 for significant difference found
                         %  0 otherwise
                    end ;

                  else ;

                    if varleft == 0 ;    %  do 1 sample t-test
                      h = ttest(dataright,meanleft,salpha) ;
                         %  1 for significant difference found
                         %  0 otherwise

                    else ;    %  do full two sample t-test

                    dfleft = length(dataleft) - 1;
										dfright = length(dataright) - 1;
										dfe  = dfleft + dfright;
										msleft = dfleft * varleft ;
										msright = dfright * varright ;
										
										difference = meanleft - meanright;
										pooleds    = sqrt((msleft + msright) * (1/(dfleft + 1) + 1/(dfright + 1)) / dfe);
										
										ratio = difference / pooleds ;
										p = 1 - tcdf(ratio,dfe) ;
										    %  one tail p-value
										p = 2 * min(p, 1-p) ;
										    %  two tail p-value

 								    h = 0;
	  								if p <= salpha, 
                        h = 1; 

                        end ;
                          %  1 for significant difference found
                          %  0 otherwise
                    
                    end ;
                
                  end ;
                    
                end ;
                  
              
                    
              elseif itest == 2 ;    %  Kolmogorov-Smirnov

                h = kstest2(dataright,dataleft,salpha) ;
                    %  1 for significant difference found
                    %  0 otherwise


              elseif itest == 3 ;    %  variances and F-test
              
                varleft = var(dataleft) ;
                varright = var(dataright) ;
            
                if varright == 0 ;
                  if varleft == 0 ;
                    h = 0 ;
                       %  no significant difference

                  else ;    %  do 1 sample t-test
                    h = 1 ;
                       %  significant difference

                  end ;

                else ;

                  if varleft == 0 ;    %  do 1 sample t-test
                    h = 1 ;
                       %  significant difference

                  else ;    %  do variance F-test

                    dfleft = length(dataleft) - 1;
										dfright = length(dataright) - 1;
	
										Fstat = varleft / varright ;
	
										p = 1 - fcdf(Fstat,dfleft,dfright) ;
										    %  one tail p-value
										p = 2 * min(p, 1-p) ;
										    %  two tail p-value

								    h = 0;
    								if p <= salpha, 
                      h = 1; 
                    end ;
                                %  1 for significant difference found
                                %  0 otherwise
                      
                  end ;
              
                end ;
              
         
              elseif itest == 4 ;    %  mean removed Kolmogorov-Smirnov

                dataleftmr = dataleft - mean(dataleft) ;
                datarightmr = dataright - mean(dataright) ;
                
                h = kstest2(datarightmr,dataleftmr,salpha) ;
                    %  1 for significant difference found
                    %  0 otherwise
                    %  Note:  have not adjusted significance for removal of mean

                  
              elseif itest == 5 ;    %  Standardized Kolmogorov-Smirnov

                dataleftst = dataleft - mean(dataleft) ;
                datarightst = dataright - mean(dataright) ;

                varleft = var(dataleftst) ;
                varright = var(datarightst) ;

                if varleft > 0 ;
                  dataleftst = dataleftst / varleft ;
                end ;

                if varright > 0 ;
                  datarightst = datarightst / varright ;
                end ;
                
                h = kstest2(datarightst,dataleftst,salpha) ;
                    %  1 for significant difference found
                    %  0 otherwise
                    %  Note:  have not adjusted significance for standardization

                    
              elseif itest == 6 ;    %  median removed Kolmogorov-Smirnov

                dataleftmr = dataleft - median(dataleft) ;
                datarightmr = dataright - median(dataright) ;
                countleft    =  length(dataleft);
                countright     =  length(dataright);
                ncoefficient      =  sqrt(2)*...
                      countleft * countright /(countleft + countright);
                
                % The sqrt(2) comes from the fact that I have 2 Brownian
                % bridges on the interval of length .5. The same applies
                % for the sqrt in the calculation of p-value.

                %              
                % Calculate F1(x) and F2(x), the empirical (i.e., sample) CDFs.
								%
								
								binEdges    =  [-inf ; sort([dataleftmr;datarightmr]) ; inf];
								
								binCounts1  =  histc (dataleftmr , binEdges);
								binCounts2  =  histc (datarightmr , binEdges);
								
								sumCounts1  =  cumsum(binCounts1)./sum(binCounts1);
								sumCounts2  =  cumsum(binCounts2)./sum(binCounts2);
								
								sampleCDF1  =  sumCounts1(1:end-1);
								sampleCDF2  =  sumCounts2(1:end-1);
								
								
								%
								% Compute the test statistic of interest.
								%
								
                deltaCDF  =  abs(sampleCDF1 - sampleCDF2);
                KSstatistic   =  max(deltaCDF);
                lambda =  (sqrt(ncoefficient) + 0.12 + 0.11/sqrt(ncoefficient))...
                                              * KSstatistic ;
                
                
                %
								%  Use the asymptotic Q-function to approximate the 2-sided P-value.
								%
                j       =  [1:101]';
                p  =  sqrt(2 * sum((-1).^(j-1).*exp(-2*lambda*lambda*j.^2)));

                h = 0;
								if p <= salpha, 
                  h = 1; 
                end ;
                            %  1 for significant difference found
                            %  0 otherwise
                            %  Note:  HAVE adjusted significance for removal of median

                  
              elseif itest == 7 ;    %  median and MAD st'dized Kolmogorov-Smirnov

                dataleftst = dataleft - median(dataleft) ;
                datarightst = dataright - median(dataright) ;

                gpmadleft = gpmad(dataleftst) ;
                gpmadright = gpmad(datarightst) ;

                if gpmadleft > 0 ;
                  dataleftst = dataleftst / gpmadleft ;
                end ;

                if gpmadright > 0 ;
                  datarightst = datarightst / gpmadright ;
                end ;
                
                h = kstest2(datarightst,dataleftst,salpha) ;
                    %  1 for significant difference found
                    %  0 otherwise
                    %  Note:  have not adjusted significance for standardization

                    
              elseif itest == 8 ;    %  Quartile st'dized Kolmogorov-Smirnov

                medianleft = median(dataleft) ;
                medianright = median(dataright) ;
                
                dataleftst = dataleft - medianleft ;
                datarightst = dataright - medianright ;

                Qoneleft =  prctile(dataleftst,25) ;
                Qoneright =  prctile(datarightst,25) ;

                Qthreeleft =  prctile(dataleftst,75) ;
                Qthreeright =  prctile(datarightst,75) ;

                flagoneleft = dataleftst < 0 ;
                flagoneright = datarightst < 0 ;

                if Qoneleft < 0 ;    %  then rescale guys on Q1 side
                  dataleftst(flagoneleft) = phiinv(0.25) * ...
                                      dataleftst(flagoneleft) / Qoneleft ;
                end ;

                if Qoneright < 0 ;    %  then rescale guys on Q1 side
                  datarightst(flagoneright) = phiinv(0.25) * ...
                                      datarightst(flagoneright) / Qoneright ;
                end ;

                if Qthreeleft > 0 ;    %  then rescale guys on Q3 side
                  dataleftst(~flagoneleft) = phiinv(0.75) * ...
                                      dataleftst(~flagoneleft) / Qthreeleft ;
                end ;

                if Qthreeright > 0 ;    %  then rescale guys on Q3 side
                  datarightst(~flagoneright) = phiinv(0.75) * ...
                                      datarightst(~flagoneright) / Qthreeright ;
                end ;
                
                countleft    =  length(dataleftst);
                countright     =  length(datarightst);
                ncoefficient      =  sqrt(4)*...
                       countleft * countright /(countleft + countright);

                % The sqrt(4) comes from the fact that I have 4 Brownian
                % bridges on the subintervals of length .25. The same applies
                % for the sqrt in the calculation of p-value.
      
                       
                %
								% Calculate F1(x) and F2(x), the empirical (i.e., sample) CDFs.
								%
								
								binEdges    =  [-inf ; sort([dataleftst;datarightst]) ; inf];
								
								binCounts1  =  histc (dataleftst , binEdges);
								binCounts2  =  histc (datarightst , binEdges);
								
								sumCounts1  =  cumsum(binCounts1)./sum(binCounts1);
								sumCounts2  =  cumsum(binCounts2)./sum(binCounts2);
								
								sampleCDF1  =  sumCounts1(1:end-1);
								sampleCDF2  =  sumCounts2(1:end-1);
								
								
								%
								% Compute the test statistic of interest.
								%
								
                deltaCDF  =  abs(sampleCDF1 - sampleCDF2);
                KSstatistic   =  max(deltaCDF);
                lambda =  (sqrt(ncoefficient) + 0.12 + 0.11/sqrt(ncoefficient))...
                                              * KSstatistic ;
                
          
                %
								%  Use the asymptotic Q-function to approximate the 2-sided P-value.
								%
                j       =  [1:101]';
                p  =  (2 * sum((-1).^(j-1).*exp(-2*lambda*lambda*j.^2)))^(0.25) ;

                h = 0;
								if p <= salpha, 
                  h = 1; 
                end ;
                            %  1 for significant difference found
                            %  0 otherwise
                            %  Note:  HAVE adjusted significance for standardization

                            
              end ;    %  end of itest if-block
              
              
              
              mout(iwid,iloc) = h + 1 ;
                    %  2 for significant difference found
                    %  1 otherwise


            end ;    %  of iloc loop

            
            
          end ;    %  of iwd loop


	
          %  Make a plot
          %
          mout = mout + 1 ;
                %  add one, to give entry numbers in colormap

          image([minloc,maxloc],[minwid,maxwid],mout) ;
            set(gca,'YDir','normal') ;
              titstr = ['ScaSta, ' dattit ', ' simultit ', ' noisetit ...
                                ', gap = ' num2str(gap) ...
                                ', ' testtit] ;
            title(titstr) ;
            xlabel('loc') ;
            ylabel('log10(wid)') ;
              
				
					if ipart == 1 ;    %  then save as a postscript file
              
              orient landscape ;
              outstr = ['\MatlabR12\work\ScaSta\outputs\ScaStaEg2' ...
                               'p' num2str(ipart) ...
                               's' num2str(isimul) ...
                               'd' num2str(idat) ...
                               'n' num2str(inoise) ...
                               't' num2str(itest) ...
                               'g' num2str(igap) ...
                                                        ] ;

            eval(['print -dpsc ' outstr]) ;  
            disp('    Finished Plot') ;
      			%gppause ;
	

    			elseif ipart == 2 ;    %  then save movie frame
      			
      			moviestruct(igap) = getframe(figh) ;

				
					end ;    %  of ipart if-block

					
					
        end ;    %  of igap loop

        
	
  			if ipart == 2 ;    %  then save and run movie

  			
  			  %  first reorder movie to make play back and forth
  			  %
          vmorder = [(1:ngap),((ngap-1):-1:2)] ;
          moviestruct = moviestruct(vmorder) ;
          %  reorder frames, to loop back to beginning
  			  
          
          %  play movie on screen
          %
				  movie(figh,moviestruct,nrepeat,fps) ;
				  
				  
				  %  save movie as avi file
				  %
          outstr = ['\MatlabR12\work\ScaSta\outputs\ScaStaEg2' ...
                           'p' num2str(ipart) ...
                           's' num2str(isimul) ...
                           'd' num2str(idat) ...
                           'n' num2str(inoise) ...
                           't' num2str(itest) ...
                                                        ] ;
				  movie2avi(moviestruct,outstr,'compression','MSVC', ...
				                               'keyframe',fps,'fps',fps) ;
  			

					
				end ;


				
      end ;    %  of inoise loop

	

    end ;    %  of idat loop

	
	
	end ;    %  of itest loop

	
	
end ;    %  of isimul if-block
