function [xn,pmean,pstd,pmci,psci]=permean1(x,per,alpha,missval,dmeanflg,normflg,plflg,datastr,f1,f2)
%
%	function [xn,pmean,pstd,pmci,psci]=permean1(x,per,alpha,missval,dmeanflg,normflg,plflg,datastr,f1,f2)
%
% x = time series vector (column)
% per = period 
% alpha = level of significance for conf intervals
% missval = code for missing data; use NaN to pass all
% dmeanflg = set >0 to remove the periodic mean from x
% normflg = set >0 to scale xn by pstd
% plflg = set to >0 to plot
% datastr = string name of data for plots; 2 plots are generated 
% pmean = the periodic mean
% pstd = periodic std dev.
% pmci = per x 2 confidence intervals for mean
% psci = per x 2 confidence intervals for std dev.
% xn = x de-meaned (if dmeanflg > 0) and normalized by pstd (if normflg > 0)
% f1 fig no
% f2 fig no.
if nargin < 9 & plflg
    f1=figure;
    f2=figure;
end
nx=length(x);
nper = floor(nx/per);
nxact=nper*per;
nrem=nx-nxact;
if plflg 
    disp(['found ',int2str(nper),' periods of length ',int2str(per),...
    ' with remaider of ',int2str(nrem)]);
end
% find and count missing in whole series
if isnan(missval)
    missisnan=1;
    imissx=find(isnan(x));
else
    missisnan=0;
    imissx=find(x==missval);
end
nmissx=length(imissx);

%  form the periodic mean and std dev ; note this way all the data is used
pmean1=zeros(nx,1);
pstd1=zeros(nx,1);
pmean=zeros(per,1);	
pstd=zeros(per,1);
ny=zeros(per,1);
pmci=zeros(per,2);
psci=zeros(per,2);
%%%
%%%   now find per mean, per std and replace missing values
%%%
for i=1:per
	index=[i:per:nx];				% index set for season i
	z=x(index);                 % data for season i
    if missisnan
        igood=find(~isnan(z));
        imiss=find(isnan(z));
    else
        igood=find(z~=missval);		% index set for good data of season i
        imiss=find(z==missval);		% missing data for season i
    end
   z=z(igood);						% good data for season i
   ny(i)=length(z);				% no. good points for season i
	pmean(i)=mean(z);				% per. mean at season i
   pstd(i)=std(z);				% per. std at season i
   x(index(imiss))=pmean(i);	% replace missing with the pmean(i)
   % form the extended periodic mean and std
   pmean1(index)=pmean(i);		% extend the mean for subtracting later
   pstd1(index)=pstd(i);			% this extends to all of x
   %  confidence intervals use statistics toolbox
  
   t0 = tinv([alpha/2 1-alpha/2],ny(i)-1);
   pmci(i,1) = pmean(i) + t0(1)*pstd(i)/sqrt(ny(i)-1);
   pmci(i,2) = pmean(i) + t0(2)*pstd(i)/sqrt(ny(i)-1);
   chi20 = chi2inv([alpha/2 1-alpha/2],ny(i)-1);
   
   psci(i,1) =  pstd(i)*sqrt((ny(i)-1)./chi20(2));
   psci(i,2) =  pstd(i)*sqrt((ny(i)-1)./chi20(1));           
   
end

xdev=x-pmean1;
if plflg
	figure(f1);
	subplot(211);
	plot(x);
	title([datastr,' original series of ',int2str(nx),' points; found ',int2str(nmissx),' missing']);
	hold on;
	plot(imissx,x(imissx),'rx');
	hold off;

	subplot(212);

	% in errorbar, the upper and lower are relative to pmean
	errorbar([1:per],pmean,pmean-pmci(:,1),pmci(:,2)-pmean,'x');
	title(['Periodic mean, No. periods = ',...
	int2str(nper),', alpha = ',num2str(alpha)]);
	v=axis;
	axis([0 per+1 v(3) v(4)]);

% now do std deviation


	figure(f2);
	subplot(211);
	plot(xdev);
	title([datastr,' deviation from periodic mean; missing(red) set to 0']);
	hold on;
	plot(imissx,xdev(imissx),'rx');
	hold off;
	subplot(212);
	errorbar([1:per],pstd,pstd-psci(:,1),psci(:,2)-pstd,'x');
	title(['Periodic standard deviations, No. periods = ',...
	int2str(nper),', alpha = ',num2str(alpha)]);
	v=axis;
   axis([0 per+1 v(3) v(4)]);
end
% now de-mean
if dmeanflg
   xn=xdev;
else
   xn=x;
end

% now normalize by pstd

if normflg
   xn = xn ./ pstd1;
end