function [mth,ym,ystd,ymin,ymax] = statsout(y,yd,fileout,ws,yr,num)
%function statsout(y,[yd],[fileout],[ws],[yr],[num])
%
%   writes statistics of matrix "y" to an output file.
%   several types of statistics are available
%
%
%      y       =  some matrix with any number of columns
%      yd      =  vector of yearday    
%      fileout =  string representing output filename
%      ws      =  a code for type of stats:
%             1 = General
%             2 = Correlation of y's columns
%             3 = Histogram
%             4 = By Month or Week with prompts for details
%             41 = By Month without prompts for details
%             42 = By Week without prompts for details 
%             5 = standard current meter & wind statistics (incl. principal axis)
%
%
%       yr       = year needed in ws=4 & 41 &42  case
%      num       = integer code to provided for the output line in
%                   order to distinguish this record from others 
%
if nargin<3;
fileout=input(' Enter output filename: ','s');
end;
%
if nargin<4;
ws=menu(' What type of stats?','General','Correlation','Histogram','By Month','current meter');
end;
%
%*********************************************************************************
if ws==5;
%
%                       STANDARDIZED CURRENT METER OUTPUT
%
% Note: In the case of current meter records "y" should be a 4 column array representing
%       temperature, u, v, and speed (or some multiple of 4)
%
% output has no header but records are as follows
%
%  fprintf(fileout,'  mean  std    min    max  ...\n')
%
%  for each of temp, e/w, n/s, and speed
%
  [nr,nc]=size(y);
  for k=1:4:nc-3;
    for ic=k:k+3;
%
%     Find observations that aren't "NaNs"
%
      f=isnan(y(:,ic));
      ff2=find(f~=1);
%
%     skip next section if ALL the records are "NaN"
%
      if length(ff2)~=0;
        ym=mean(y(ff2,ic));
        ystd=std(y(ff2,ic));
        ymin=min(y(ff2,ic));
        ymax=max(y(ff2,ic));
%
        fprintf(fileout,'%8.2f%6.2f',ym,ystd)
%
%       In the case of speed do not report minimum 
%
        if rem(ic,4)==0;
          fprintf(fileout,'%6.1f',ymax);
        else;
          fprintf(fileout,'%7.1f%6.1f',ymin,ymax);
        end;
      end;
    end
%
%   Finally, put out the prinicpal axis angle (COUNTERclockwise from east)
%   and add an end of line code (don't plot)
%

    [ang,xs,ys]=prinaxis2(y(ff2,k+2),y(ff2,k+3),0);
    fprintf(fileout,'%7.1f %7.1f %7.1f\n',ang,xs,ys);
  end;
end
%
%***************************************************************************
if ws==1
%
%                     BASIC STATS
%

fprintf(fileout,'\n variable  mean  std    min    max ...\n')
[nr,nc]=size(y);
  for i=1:nc;
    ym=mean(y(:,i));
    ystd=std(y(:,i));
    ymin=min(y(:,i));
    ymax=max(y(:,i));
%
      fprintf(fileout,'%9.0f',i);
      fprintf(fileout,'%7.2f%7.2f',ym,ystd)
      fprintf(fileout,'%7.1f%7.1f\n',ymin,ymax)
  end;
end
%
%***************************************************************************
if ws==2
%
%                     CORRELATION
%
% write out the correlation matrix
  cc=corrcoef(y);
  save stats.out cc /ascii
end
%
%**************************************************************************
if ws==3
%
%                   HISTOGRAM 
%
% write out the histogram stats per bin
  [m,p]=size(y);
  for i=1:p
    fprintf(fileout,' Variable: %2.0f',i);
    [n,x]=hist(y(:,i));
    np=n';
    xp=x';
    h=[xp,np];
    for j=1:length(h);
     fprintf(fileout,'\n%f%f',h(j,1),h(j,2));
    end
  end
end
%
%**********************************************************************
if ws ==4|ws==41|ws==42;
%
%                   prints MEANS over some time period (weekly,mthly,etc)
%
           disp('Hardcoded to expect u & v in 2nd and 3rd columns,')
           disp('respectively, for calculation of resultant.')


  if nargin<2;
    disp(' ');
    disp(' Option requires second argument to be passed to "statsout" function,')
    disp('  a yearday vector.')
    disp(' ')
  end;
%
% "ms" is the month (or week) start yearday
%
  if ws==41;
    morw=1;
    disp('Hardcoded to generate MONTHLY stats');
  elseif ws==42;
    morw=2;
    disp('Hardcoded to generate MONTHLY stats');
  else;
    morw=menu('What?','by month','by week');
  end;

  if morw==1;
    ms=[0 31 59 90 120 151 180 212 243 273 304 334 365];
  else;
    ms=[0:floor(365/52):365];
  end;


  mstr=['jan';'feb';'mar';'apr';'may';'jun';'jul';'aug';'sep';'oct';'nov';'dec'];
%
%  fprintf(fileout,'\n   month     mean    std     min    max   prinaxis  npts...\n')
%
[nr,nc]=size(y);
%
% modify the "ms" vector for leap years
%
  if rem(yr,4)==0;
    idf=find(ms>59);
    ms(idf)=ms(idf)+1;
  end;

% decide whethter to calculate & report prinicpal axis
  if ws==4;
    pa=menu('Do you want principal axis stats?','yes','no');
  else;
    pa=2;
    disp('Hardcoded to NOT generate principal axis stats');
    disp(' ');
  end;

%
% Begin loop for 12 months (or 52 weeks)
%
  for mm=1:length(ms)-1;
    ff=find(yd>=ms(mm)&yd<ms(mm+1));
    if length(ff)~=0;
%
%   Begin loop for case where "y" is a single variable (ie. one column vector)
%
     if nc==1;
%
%    Find observations that aren't "NaNs"
%
       f=isnan(y(ff));
       ff2=find(f~=1);
%
% skip next section if ALL the records are "NaN"
%
       if length(ff2)~=0;
        mth(mm)=mm; 
        ym(mm)=mean(y(ff(ff2)));
        ystd(mm)=std(y(ff(ff2)));
        ymin(mm)=min(y(ff(ff2)));
        ymax(mm)=max(y(ff(ff2)));
        npts=length(y(ff(ff2)));
%
% write out results
%
        fprintf(fileout,'%3.0f %3.0f',yr,mm);
        fprintf(fileout,'%8.3f%7.2f',ym(mm),ystd(mm))
        fprintf(fileout,'%7.1f%7.1f%7.0f\n',ymin(mm),ymax(mm),npts);
       end;
%
    else;
%
%   Begin loop for "nc" variables
%
      for ic=1:nc;
%
%    Find observations that aren't "NaNs"
%
         f=isnan(y(ff,ic));
         ff2=find(f~=1);
%
% skip next section if ALL the records are "NaN"
%
        if length(ff2)~=0;
          mth(mm)=mm; 
          ym(mm)=mean(y(ff(ff2),ic));
%
% save this component for resultant calculation
% assumes u and v are 2nd  and 3rd variable
%
           if ic==2;
             umean=ym(mm);
           elseif ic==3;
             vmean=ym(mm);
           end;
%
        ystd(mm)=std(y(ff(ff2),ic));
        ymin(mm)=min(y(ff(ff2),ic));
        ymax(mm)=max(y(ff(ff2),ic));
        npts=length(y(ff(ff2),ic));
%
% write out results
% 
% (If an integer code is provided (nargin=6) then put it as first line)
%
        if ic==1;
          if nargin>5;
            fprintf(fileout,'%5.0f%3.0f%3.0f',num,yr,mm);
          else;
            fprintf(fileout,'%3.0f %3.0f',yr,mm);
          end;
        end;
        fprintf(fileout,'%8.3f%7.2f',ym(mm),ystd(mm))
        fprintf(fileout,'%7.1f%7.1f',ymin(mm),ymax(mm));
%
% if this is the last variable, calculate resultant spd & dir and principal axis.
% put a carriage return and # of pts only if this is the last column of the y matrix
%
        if ic==nc;


% calculation taken from "directio.m" routine and,
% if the user requested, the "prinaxis2" routine

          [dir,spd]=directio(umean,vmean);
          if pa==1;
            [deg,xs,ys]=prinaxis2(y(ff(ff2),2),y(ff(ff2),3),0);
            fprintf(fileout,'%6.0f%9.3f%5.0f',npts,spd,dir);
            fprintf(fileout,'%5.0f%8.3f%9.3f\n',deg,xs,ys);
          elseif pa==2;
            fprintf(fileout,'%6.0f %9.3f %5.0f\n',npts,spd,dir);
          end;
        end;

% if length(ff2~=0)
        end;

% for ic=1:nc;
       end;

% if nc==1;
      end;

% if length(ff)~=0
    end;

% for m=1:12;
  end;

% if ws==4; 
end;  

