% this script makes a batch of par simulations
% the role of makepar1 is included and makepar is called
% in a loop of nr simulations
% parameters are n,T,ph0,phidel,phiph,del0,deldel,delph
%par(1,:)=[512*50 32 .1 .5 1.5*pi 2 0 0];
%par(2,:)=[1000 12 .6 .3 1.5*pi 2 0 0];
%par(3,:)=[1000 12 .9 .3 1.5*pi 2 0 0];
%PRINTSET=[2 3 4 5 9 10];
PRINTSET=[];
numpts=512*10;
par(1,:)=[numpts 32 .2 .2 1.5*pi 1 0 0];
par(2,:)=[numpts 32 .4 .2 1.5*pi 1 0 0];
par(3,:)=[numpts 32 .8 .2 1.5*pi 1 0 0];
par(4,:)=[numpts 32 .4 .4 1.5*pi 1 0 0];
par(5,:)=[numpts 32 .8 .4 1.5*pi 1 0 0];
par(6,:)=[numpts 32 .8 .6 1.5*pi 1 0 0];
par(7,:)=[numpts 32 .8 .8 1.5*pi 1 0 0];
par(8,:)=[numpts 32 -.2 .2 1.5*pi 1 0 0];
par(9,:)=[numpts 32 -.4 .2 1.5*pi 1 0 0];
par(10,:)=[numpts 32 -.8 .2 1.5*pi 1 0 0];
par(11,:)=[numpts 32 -.4 .4 1.5*pi 1 0 0];
par(12,:)=[numpts 32 -.8 .4 1.5*pi 1 0 0];
par(13,:)=[numpts 32 -.8 .6 1.5*pi 1 0 0];
par(14,:)=[numpts 32 -.8 .8 1.5*pi 1 0 0];
% the following were simulations based on parameters estimated
% when par1 had an error
%par(1,:)=[8112 24 .277 .058 12.6*pi/180 3.89 .28 -pi];
%par(2,:)=[8112 24 .099 .06 137*pi/180 3.35 .19 125*pi/180];
%par(3,:)=[8112 24 .358 .096 63*pi/180 2.13 .12 76*pi/180];
%par(4,:)=[8112 24 .373 .044 -80*pi/180 3.35 .29 27*pi/180];
%par(5,:)=[8112 24 .139 .092 117*pi/180 4.42 1.14 147*pi/180];
% the following aare based on the parameters estimated
% from the txozone data using the corrected par1
%par(1,:)=[8112 24 .9 .035 52*pi/180 1.73 .17 50*pi/180];
%par(2,:)=[8112 24 .82 .027 -129*pi/180 1.94 .15 101*pi/180];
%par(3,:)=[8112 24 .84 .032 -27*pi/180 1.36 .14 -168*pi/180];
%par(4,:)=[8112 24 .89 .04 -35*pi/180 1.62 .32 122*pi/180];
%par(5,:)=[8112 24 .85 .085 -153*pi/180 2.33 .33 -70*pi/180];
[nr,nc]=size(par);
for k=1:nr
	n=par(k,1);
	T=par(k,2);
	ph0=par(k,3);
	phidel=par(k,4);
	phiph=par(k,5);
	del0=par(k,6);
	deldel=par(k,7);
	delph=par(k,8);
	prd=1;
	for i=1:T
		phi(i,1)=ph0+phidel*cos((2*pi*i/T)+phiph);
		del(i,1)=del0+deldel*cos((2*pi*i/T)+delph);
		prd=prd*phi(i,1);
   end
   basename=['parsim',sprintf('%0.2d',k)];
	figure(1);
	subplot(211);
	plot(phi(:,1));
	mod=phidel/ph0;
   title(['phi  - deviation to mean: ',num2str(mod),' for ',basename]);
   subplot(212);
	plot(del(:,1));
	mod=deldel/del0;
   title(['del  - deviation to mean: ',num2str(mod),' for ',basename]);
   
	disp([basename,'  ph0:',num2str(ph0),' phidel:',num2str(phidel),' product phi(t): ',...
	num2str(prd)]); 
	x(k,:)=makepar(n,phi,del);

	[xn,pmean,pstd,pmci,psci]=permean1(x(k,:)',T,.05,NaN,1,0,0,basename);
   %  dmean=1,normflg=0, plflg=0
   figure(2);
   subplot(211);
	plot(x(k,1:256));
   title(['first 256 points of case: ',basename]);
   v=axis;
   axis([1 256 v(3) v(4)]);
   
   subplot(212);
   cv=par1cov(phi',del',1);
   errorbar([1:T],pstd,pstd-psci(:,1),psci(:,2)-pstd,'x');
   hold on;
   plot(sqrt(diag(cv)));
   title(['theoretical and estimated standard deviation - one period']);
   v=axis;
   axis([1 T v(3) v(4)]);
   hold off;

   % sqplot will generate fig 3
   sqplot(x(k,:)',2,80,512,1,1,.01,8,basename);
   
  
	
   if ismember(k,PRINTSET)
      for j=1:3
         figure(j);
         switch j
         case 1
            name=[basename,'pm'];
         case 2
            name=[basename,'ts'];
         case 3
            name=[basename,'sq'];
         end
         print('-deps',['.\fourteen\',name]); % put in subdir called fourteen
      end
   end
	pause;
	close all;
end
%fout= fopen(['c:\mstuff\par\pbat001.s'],'w');

% fwrite goes in column order
%fwrite(fout,y','float');		

%fclose(fout);