%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
	out1=sprintf('\n\n wavelet_synth_alt.m,  ver 1.0, April 8, 2011 \n\n by Tom Irvine \n Email: tomirvine@aol.com \n');
	out2=sprintf(' This program synthesizes a time history using wavelets to satisfy ');
	out3=sprintf(' a shock response spectrum (SRS) specification. \n');
%
	out4=sprintf('\n The program also optimizes the time history to yield the lowest overall \n error, acceleration, velocity, and displacement. \n');
    out5=sprintf('\n The optimization is performed via trial-and-error.\n');
%
    disp(out1);
    disp(out2);
    disp(out3);
    disp(out4);
    disp(out5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
MAXTRIALS=20000;
K100=100000;
NUM=125;
MAX=125;
exponent=1;  % initialize only
%
tpi=2.*pi;
%
niter=100;
%
clear nspec;
clear store_amp;
clear store_NHS;
clear store_td;
%
clear nhs;
%
clear ym;
clear vm;
clear dm;
clear em;
clear im;
clear sym;
clear svm;
clear sdm;
clear sem;
clear sim;
%
clear yrank;
clear vrank;
clear drank;
clear erank;
clear irank;
%
clear pyrank;
clear pvrank;
clear pdrank;
clear perank;
clear pirank;
%
clear nrank;
clear upper;
clear a1;
clear a2;
clear b1;
clear b2;
clear b3;
%
clear local_amp; 
clear amp;
%
clear f;
clear fn;
clear fr;
%
clear omegaf[MAX];
%
clear over_period;
clear onep5_period;
%
clear r;
clear spec;
clear s;
clear amp_start;
clear td;
clear th;
%
clear tmax;
clear tmin;
%
clear vth;
clear wavelet;
%
clear x;
clear xb;
clear xbb;
clear xmax;
clear xmin;
clear xxmin;
clear xxmax;
%
clear yb;
clear ybb;
%
%%	out6=sprintf('\n Select SRS algorithm:');
%%	out7=sprintf('\n 1=Kelly-Richman  2=Smallwood \n');
%
%%    disp(out6);
%%    disp(out7);
%
%%    algrorithm=input(' ');
%
    algorithm=1; 
%
     th=zeros(K100,1);
    vth=zeros(K100,1);    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    clear f;
    clear a;
%
    disp(' Select data input method.');
    disp('  1=keyboard  ');
    disp('  2=internal Matlab array  ');
    disp('  3=external ASCII file   ');
	imethod=input(' ');
    disp(' ');
%	  
    if(imethod==1)
%
        n=0;
        while(n < 2)
            disp('  How Many Breakpoints (min=2)? ');
            n=input(' ');
        end
%			  
        disp('  Enter Points (free-format)');
        disp(' ');
%       
        f=zeros(n,1);
        a=zeros(n,1);
        for i=1:n
                out1=sprintf(' Enter natural frequency(Hz) %d:  ',i);
                disp(out1);
                f(i)=input(' ');
                out2=sprintf(' Enter SRS(G) %d:  ',i);
                disp(out2);
                a(i)=input(' ');  
                disp(' ')
        end
%	
    end
%
    if(imethod==2)
        disp('The array must have two columns:  Natural Freq(Hz)  SRS(G)   ');
        THM=input(' Enter the array name:  ');
        f=THM(:,1);
        a=THM(:,2);
        clear length;
        n=length(a);
    end
%
    if(imethod==3)
%
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename); 
        fid = fopen(filename,'r');
        THF = fscanf(fid,'%g %g',[2 inf]);
        THF=THF';   
        f=THF(:,1);
        a=THF(:,2);
        clear sz;
        sz=(size(f));
        n=sz(1);
        if(sz(2)>sz(1))
            n=sz(2);
        end
%
    end
    clear aspec;
    aspec=a;
%
    ffirst=f(1);
    flast=f(n);
    last_f=f(n);
    last_a=a(n);
%
    clear fr;
    clear r;
    fr=f;
    r=a;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%  Calculate slopes between input points
%
      clear length;
      num=length(fr);
%
      s=zeros(num-1,1);
%
       for i=1:(num-1)
          a=(log(r(i+1))-log(r(i)));
          b=(log(fr(i+1))-log(fr(i)));
          s(i)=a/b;
       end
%
   clear f;
   clear spec;
%
%%  Interpolate
%
	   f(1)= fr(1);
	   fa=fr(1);
       fb=fr(1);
	   spec(1)=r(1);
	   amp_start(1)=(spec(1)/16.);
%       
	   i=1;
%
	   out1=sprintf('\n\n Enter octave spacing. \n 1= 1/3   2= 1/6   3= 1/12 \n');
       disp(out1);
%	   
       ioct=input(' ');
%
	   octave=(1./3.);
%
	   if ioct==2
		  octave=(1./6.);
       end
%
	   if ioct==3
		  octave=(1./12.);
       end
%
       f(1)=fr(1);
       spec(1)=r(1);
       i=2;
	   while(1) 
		  ff=(2.^octave)*fb;
		  fb=ff;
		  if(ff>fr(num))
              break;
          end    
%
		  if( ff >= fr(1))
%
				for j=1:num
%
					if(ff == fr(j))
						f(i)=ff;
						spec(i)=r(j);
                        nspec=i;
                        i=i+1;                        
						break;
                    end
					if(ff < fr(j) && j>1)
%					
						f(i)=ff;
						az=(log10(r(j-1)));
						az=az+(s(j-1)*(log10(ff)-log10(fr(j-1))));
						spec(i)=10.^az;
                        nspec=i;
                        i=i+1;
						break;
                    end
%
                end
          end
       end
%
       frlast=max(fr);
%
	   if(frlast > f(nspec))
          nspec=nspec+1;
		     f(nspec)=fr(num);
		  spec(nspec)=r(num);
       end
	   amp_start= spec/16;
       arlast=amp_start(nspec);   
%
	   if(nspec > NUM)  
		   out1=sprintf('\n Warning: number of specification points reduced. ');
           disp(out1);
		   nspec=NUM;
       end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
	  out1=sprintf('\n Enter damping format for SRS. \n 1= damping ratio   2= Q \n');
      disp(out1);
%
	  idamp=input(' ');
%
      while(1)
%
        if(idamp == 1 )
            disp(' Enter SRS damping ratio (typically 0.05) ');
            damp=input(' ');
        else
            disp(' Enter SRS amplification factor Q (typically 10) ');
            damp=input(' ');
            damp = (1./(2.*damp));
        end
        if(damp<1.0)
            break;
        end
%
      end
%
      disp(' ');
	  disp(' Enter the number of trials. ');
%      
      ntrials=input(' ');
%
	  if(ntrials> MAXTRIALS)
%	  
		  ntrials=MAXTRIALS;
%
          out1=sprintf('\n Warning: number of trials reduced to %ld \n',ntrials);
          disp(out1);
      end
%
    ym=zeros(ntrials,1);
    vm=zeros(ntrials,1);
    dm=zeros(ntrials,1);
    em=zeros(ntrials,1);   
    im=zeros(ntrials,1);      
%
	  while(1)
%         
		out1=sprintf('\n\n Enter units ');  
		out2=sprintf('    1=English:  G,       in/sec, in   ');
		out3=sprintf('    2=metric:   G,       m/sec,  mm   ');
		out4=sprintf('    3=metric:   m/sec^2, m/sec,  mm  \n');
%
        disp(out1);
        disp(out2);
        disp(out3);
        disp(out4);
%
        iunit=input(' ');
%        
        if( iunit==1 || iunit==2 || iunit==3 )
            break;
        end
      end
%
	  dunit='mm';
	  vunit='m/sec';
	  aunit='G';
%
	  if(iunit==1)
		  dunit='inch';
		  vunit='in/sec';
      end
	  if(iunit==3)
		  aunit='m/sec^2';
      end
%
      if(iunit==1)
          disp(' ');
          disp(' Enter displacement limit (in zero-to-peak) ');
          displacement_limit=input(' ');
      else
          disp(' ');
          disp(' Enter displacement limit (mm zero-to-peak) ');
          displacement_limit=input(' ');       
      end
%
      if(iunit==1)
          disp(' ');
          disp(' Enter velocity limit (in/sec) ');
          velocity_limit=input(' ');
      else
          disp(' ');
          disp(' Enter velocity limit (m/sec) ');
          velocity_limit=input(' ');       
      end    
%
      if(iunit==3)
          disp(' ');
          disp(' Enter acceleration limit (m/sec^2) ');
          acceleration_limit=input(' ');
      else
          disp(' ');
          disp(' Enter acceleration limit (G) ');
          acceleration_limit=input(' ');       
      end        
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    frmax=max(fr);
%
	sr=(10*frmax);
%
	out1=sprintf(' Recommend sample rate >= %10.6g samples/sec',sr );
    out2=sprintf('\n Enter sample rate ');
%
    disp(out1);
    disp(out2);
%
	sr=input(' ');
%
	if(sr < 4*frmax)
		out1=sprintf('\n Warning: sample rate reset to %f ',sr);
        disp(out1);
		sr= 4*frmax;
    end
%
	dt=(1./sr);
	dth=(dt/2.);
%
	out1=sprintf('\n\n Select condition  ');
	out2=sprintf('   1=limit time history points   2=set duration \n');
%
    disp(out1);
    disp(out2);
%
    idur=input(' ');
%
	if(idur==2)
%        
		out1=sprintf('\n\n Enter duration (sec) ');
		out2=sprintf(' (Recommend %6.3f or greater) \n',2.0/f(1));
%
        disp(out1);
        disp(out2);
%
        dur=input(' ');
%
	    nt=round(dur/dt);        
%
		if(nt>K100)
			out1=printf('\n\n Warning: duration reduced. \n');
            disp(out1);
			nt=K100;
        end 
%
    else
%
		out1=sprintf('\n Enter the number of points \n');
		disp(out1);
        nt=input(' ');
%
		if(nt>K100)
			out1=sprintf('\n\n Warning: duration reduced. \n');
            disp(out1);
			nt=K100;
        end
%
		dur=nt*dt;
%
    end
%
	out1=sprintf('\n dt=%9.4g sec   dur=%8.4f sec  sr=%9.4g sample/sec  nt=%ld \n',dt,dur,sr,nt);
    disp(out1);
%
	if(dur < 1.5/f(1))
%
        dur=1.5/f(1);
		out1=sprintf('\n\n Warning: duration is too short.\n\n Duration is reset to %f ',dur);
        disp(out1);
%
    end
%        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' The following weight numbers will be used to select the optimum waveform. ');
disp(' Suggest using integers from 0 to 10 ');
disp(' ');
iw=input(' Enter individual error weight ');
ew=input(' Enter      total error weight ');
dw=input(' Enter     displacement weight ');
vw=input(' Enter         velocity weight ');
aw=input(' Enter     acceleration weight ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
	disp(' ');
%
    a1=zeros(nspec,1);
    a2=zeros(nspec,1);
    b1=zeros(nspec,1);
    b2=zeros(nspec,1);   
    b3=zeros(nspec,1);     
%
    E=0;
    K=0;
    C=0;
    S=0;
    Sp=0;
%
	for j=1:nspec
%
		omega=(tpi*f(j));
		omegad=(omega*sqrt(1.-damp^2));
%            
		if(algorithm==1)
%
			cosd=(cos(omegad*dt));
			sind=(sin(omegad*dt));
			domegadt=(damp*omega*dt);

			a1(j)=(2.*exp(-domegadt)*cosd);
			a2(j)=(-exp(-2.*domegadt));
			b1(j)=(2.*domegadt);
			b2(j)=(omega*dt*exp(-domegadt));
			b2(j)=b2(j)*(( (omega/omegad)*(1.-2.*(damp^2))*sind -2.*damp*cosd ));
			b3(j)=0.;
%           
        else
%
			E=(exp(-damp*omega*dt));
			K=(omegad*dt);
			C=(E*cos(K));
			S=(E*sin(K));
%
			Sp=S/K;
%
			a1(j)=(2*C);
			a2(j)=(-(E^2));
%
			b1(j)=(1.-Sp);
			b2(j)=(2.*(Sp-C));
			b3(j)=((E^2)-Sp);
        end
    end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
		errorbefore=0.;
%
		ijk=0;
%
		if(nt < K100)
%
			    out1=sprintf('             Peak        Peak         Peak         ');
				out2=sprintf('   Trial     Accel       Vel          Disp      T.Error      I.Error');
%
                disp(out1);
                disp(out2);
%                
				if(iunit==1)
					          out3=sprintf('             (G)         (in/sec)     (in)       ');
                end
				if(iunit==2)
					          out3=sprintf('             (G)         (m/sec)      (mm)       ');
                end
				if(iunit==3)
					          out3=sprintf('             (m/sec^2)   (m/sec)      (mm)       ');
                end
                disp(out3);
%
            ichoice = 1;
%
            omegaf=tpi*f;
            local_amp=zeros(nspec,1);
            over_period=zeros(nspec,1);
            onep5_perios=zeros(nspec,1);
%
			for i=1:nspec 
				over_period(i)=(1.0/f(i));
			   onep5_period(i)=(1.5/f(i));
            end
%
			limit=round(1.2/(dt*f(1)));
%
            wavelet=zeros(nspec,nt);
%
			for inn=1:ntrials			
				local_record = 1.0e+99;
%
				exponent = 0.5;
%
				if( rand()< 0.4 )
%
					exponent = ( 0.40 + 0.1*rand() );
%
                end
%
%%				if(rand()<0.5 )
%%                    [amp,td,nhs,stype]=ws_synth3(amp_start,dur,onep5_period,inn,ntrials,f,nspec);
%%               else
%%                  [amp,td,nhs,stype]=ws_synth1(amp_start,dur,onep5_period,inn,f,nspec);
%%                end                    
                [amp,td,nhs,stype]=ws_synth3(amp_start,dur,onep5_period,ntrials,inn,f,nspec);
                 nspec=max(size(spec));
%
				ichoice=ichoice+1;
				igen=inn;
%
				for i=1:nspec
					upper(i) = (nhs(igen,i)/(2.*f(i)));
                    if( nhs(igen,i) == 0 || abs(upper(i))<1.0e-20  )
                         out1=sprintf(' f=%8.4g  u=%8.4g  nhs=%d',f(i),upper(i),nhs(igen,i));
                         disp(out1);
                        input(' ctrl-C');
                    end
                end  
%%%
%
                        re1=0;
                        re2=0;
                        rnv=0;
				for nv=1:niter  
%
				    [wavelet,th] = ws_gen_time(nhs,amp,omegaf,upper,nt,dt,td,igen,nv,wavelet,nspec);  
%                    
					[xxmax,xxmin,xmax,xmin]=ws_srs(nspec,th,a1,a2,b1,b2,b3,f);
%                    
					[error,irror]=ws_srs_error(spec,xmax,xmin,nspec);
%
					if( irror < local_record )   %% end criteria
%					
						for(i=1:nspec)
%						
							local_amp(i)=amp(igen,i);
							local_record = irror;
%
                        end
                    end
%				
%%                    out1=sprintf('%d  %8.4g',nv,irror);
%%                    disp(out1);
					if(nv>= 40 && irror >= errorbefore)
                        break;
                    end
%
                    if(nv>= 2 && irror < errorbefore)
                        %
						for i=1:nspec
							amp(igen,i)=local_amp(i);
        		        end
%
       		            store_amp(igen,:)=amp(igen,:);
      		            store_NHS(igen,:)=nhs(igen,:);
      		            store_td(igen,:) =td(igen,:);
%
     		            ymax=max(abs(th));
						[vmax,dmax]=ws_max_param(iunit,th,nt,dt);
%
						irror=(local_record);
%
                        re1=20.*error;
                        re2=20.*irror;
                        rnv=nv;
%
					    sym(inn)=abs(ymax);
					    svm(inn)=abs(vmax);
						sdm(inn)=abs(dmax);
						sem(inn)=abs(error);
						sim(inn)=abs(irror);
                    end
%
					errorbefore=irror;
%
					[amp]=ws_scale(xmax,xmin,spec,exponent,amp,inn,nspec);
%
                end
                out1=sprintf('      %ld %10.2f %12.3f %12.3f %11.2f %10.2f  %ld ',inn,ymax,vmax,dmax,re1,re2,rnv);
     		    disp(out1);
%
				ijk=ijk+1;
            end
%
			for jk=1:ntrials
				ym(jk)=sym(jk);
				vm(jk)=svm(jk);
				dm(jk)=sdm(jk);
				em(jk)=sem(jk);
				im(jk)=sim(jk);
            end
%
			rntrials=ntrials;
%
           [iwin,nrank,dm,drank]=ws_rankfunctions_alt(rntrials,ym,vm,dm,em,im,nspec,aunit,vunit,dunit,displacement_limit,velocity_limit,acceleration_limit,iw,ew,dw,vw,aw);
%
           if(iwin==0)
               disp(' Error:  no case found.  Try increasing displacement limit. ');
           else
%
           [acceleration,velocity,displacement]=ws_th_from_wavelet_table(iwin,store_amp,store_NHS,store_td,dur,dt,nt,iunit,f);
%
           figure(1);
           plot(acceleration(:,1),acceleration(:,2));
           grid on;
           title('Acceleration');
           xlabel('Time (sec)');
           if(iunit==3)
              ylabel('Accel (m/sec^2)');
           else
              ylabel('Accel (G)');               
           end
                             h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11)
%          
%
           figure(2);
           plot(velocity(:,1),velocity(:,2));
           xlabel('Time (sec)');
           if(iunit==2 || iunit==3)
              ylabel('Velocity (m/sec)');
           else
              ylabel('Velocity (in/sec)');               
           end
           grid on;
           title('Velocity');
                             h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11)
%          
%
           figure(3);
           plot(displacement(:,1),displacement(:,2));
           xlabel('Time (sec)');
           if(iunit==2 || iunit==3)
              ylabel('Disp (mm)');
           else
              ylabel('Disp (in)');               
           end
           grid on;
           title('Displacement');
                             h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11)
%          
           clear wavelet_table;
           for(i=1:nspec)
             wavelet_table(i,1)=i;
             wavelet_table(i,2)=store_amp(iwin,i);
             wavelet_table(i,3)=f(i);
             wavelet_table(i,4)=store_NHS(iwin,i);
             wavelet_table(i,5)=store_td(iwin,i);
           end
%          
           disp(' Output Time Histories:');
           disp('   displacement ');
           disp('   velocity ');
           disp('   acceleration ');
           disp('   shock_response_spectrum ');
           disp(' ');
           disp('   wavelet_table  [index  accel(G)  freq(Hz)  number of half-sines  delay(sec)]');
%           
           th=acceleration(:,2);
           [xxmax,xxmin,xmax,xmin]=ws_srs(nspec,th,a1,a2,b1,b2,b3,f);           
%
           figure(4);
           plot(f,xmax,'blue',f,xmin,'red',fr,aspec,'black');
           clear Shock_Response_Spectrum;
           Shock_Response_Spectrum=[f',xmax,xmin];
           if(iunit==3)
              ylabel('Peak Accel (m/sec^2)');
           else
              ylabel('Peak Accel (G)');               
           end
           xlabel('Natural Frequency (Hz)');
           legend ('positive','negative','spec & tol',3);  
           set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');  
           clear Q;
           Q = 1/(2*damp);
           out5 = sprintf(' Shock Response Spectrum Q=%g ',Q);
           title(out5);   
           grid on;
           hold on;
           plot(fr,aspec/sqrt(2),'black',fr,aspec*sqrt(2),'black');
           hold off;
%
%
          srs_max=max(aspec*sqrt(2));
              srs_min=min(aspec/sqrt(2));
              ymax= 10^(round(log10(srs_max)+0.8));
              ymin= 10^(round(log10(srs_min)-0.6));
%
              fmax=max(fr);
              fmin=fmax/10.;
%
              fmax= 10^(round(log10(fmax)+0.5));
%
             if  fr(1) >= 0.001
                  fmin=0.001;
              end
              if  fr(1) >= 0.01
                  fmin=0.01;
              end
              if  fr(1) >= 0.1
                  fmin=0.1;
              end
              if  fr(1) >= 1
                  fmin=1;
              end
              if  fr(1) >= 10
                  fmin=10;
              end
              if  fr(1) >= 100
                  fmin=100;
              end
%
              if(round(max(fr))==10000)
                  fmax=2000;
              end
              if(round(max(fr))==1000)
                  fmax=2000;
              end             
              if(round(max(fr))==2000)
                  fmax=2000;
                  if(round(min(fr))==10)
                          set(gca,'xtick',[10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 2000])
                          set(gca,'XTickLabel',{'10';'';'';'';'';'';'';'';'';'100';'';'';'';'';'';'';'';'';'1000';'2000';})
                  end
              end
              axis([fmin,fmax,ymin,ymax]);
                  h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11)
%          
           end
        else
			out1=sprintf('\n\n  Error: too many samples.\n\n');
            disp(out1);
        end  %% end nt loop
disp(' '); 
toc