disp('Running MATLAB script file gpfoant.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION gpfoan,
%    General Purpose Family + SiZer Approach to Fourier Analysis


idat = 0 ;     %    0 - to loop through all simple tests
               %    1 - White Noise
               %    2 - Low Frequency sin
               %    3 - High Frequency sin
               %    4 - Low Frequency sin + noise
               %    5 - High Frequency sin + noise

               %  100 - loop through variations on the sunspot data
               %  101 - Sunspot Data
               %  102 - Sunspot Data - mean
               %  103 - sqrt(Sunspot Data)
               %  104 - sqrt(Sunspot Data) - mean

               %  110 - loop through variations on the Lynx data
               %  111 - Canadian Lynx Data
               %  112 - Canadian Lynx Data - mean
               %  113 - log10(Canadian Lynx Data)
               %  114 - log10(Canadian Lynx Data) - mean


if idat == 0 ;
  idats = 1 ;
  idate = 5 ;
elseif  idat == 100 ;
  idats = 101 ;
  idate = 104 ;
elseif  idat == 110 ;
  idats = 111 ;
  idate = 114 ;
else ;
  idats = idat ;
  idate = idat ;
end ;



if idat < 100 ;    %  then work with simple tests

  for iidat = idats:idate ;

    clf ;
    if iidat == 1 ;
      etitstr = 'gpfoan test 1: White Noise' ;
      n = 256 ;
      seed = 92873489 ;
      randn('seed',seed) ;
      data = randn(n,1) ;

    elseif iidat == 2 ;
      etitstr = 'gpfoan test 2: Low Frequency sin' ;
      n = 256 ;
      data = sin(2 * pi * (1:n)' / 32) ;

    elseif iidat == 3 ;
      etitstr = 'gpfoan test 3: High Frequency sin' ;
      n = 256 ;
      data = sin(2 * pi * (1:n)' / 4) ;

    elseif iidat == 4 ;
      etitstr = 'gpfoan test 4: Low Frequency sin + Noise' ;
      n = 256 ;
      data = sin(2 * pi * (1:n)' / 32) ;
      seed = 75934753 ;
      randn('seed',seed) ;
      sig = 0.2 ;
      data = data + sig * randn(n,1) ;


    elseif iidat == 5 ;
      etitstr = 'gpfoan test 5: High Frequency sin + Noise' ;
      n = 256 ;
      data = sin(2 * pi * (1:n)' / 4) ;
      seed = 50875049 ;
      randn('seed',seed) ;
      sig = 0.2 ;
      data = data + sig * randn(n,1) ;

    end ;


    gpfoan(data,0,1,0,1) ;
          %  These are chosen as:
          %    0 - taper prop = 0, no tapering 
          %    1 - icolor = 1, full color
          %    0 - nhp, use default of 11
          %    ptrans (0 for log, 1 for raw)
      subplot(3,1,1) ;
        title(etitstr) ;


    if idat == 0 ;    %  then are looping, so pause
      disp(['Any key to continue']) ;
      pause ;
    end ;


  end ;

else ;    %  then work with real data sets

  for iidat = idats:idate ;

    clf ;
    if iidat == 101 ;
      dtitstr = 'Sunspots' ;
      infstr = 'sunspot.dat' ;
      frmtstr = '%g' ;
      insize = [1, inf] ;
      vicol = [1] ;
      datatran = [] ;
      xstr = 'Year' ;
      ystr = 'Sunspots' ;
    elseif iidat == 102 ;
      dtitstr = 'Sunspots - mean' ;
      infstr = 'sunspot.dat' ;
      frmtstr = '%g' ;
      insize = [1, inf] ;
      vicol = [1] ;
      datatran = 'data = data - mean(data) ;' ;
      xstr = 'Year' ;
      ystr = '# SS - mean' ;
    elseif iidat == 103 ;
      dtitstr = 'sqrt(Sunspots)' ;
      infstr = 'sunspot.dat' ;
      frmtstr = '%g' ;
      insize = [1, inf] ;
      vicol = [1] ;
      datatran = 'data = sqrt(data) ;' ;
      xstr = 'Year' ;
      ystr = 'sqrt(# Sunspots)' ;
    elseif iidat == 104 ;
      dtitstr = 'sqrt(Sunspots) - mean' ;
      infstr = 'sunspot.dat' ;
      frmtstr = '%g' ;
      insize = [1, inf] ;
      vicol = [1] ;
      datatran = 'data = sqrt(data) ; data = data - mean(data) ;' ;
      xstr = 'Year' ;
      ystr = 'sqrt(# SS) - mean' ;

    elseif iidat == 111 ;
      dtitstr = 'Canadian Lynx' ;
           %  Famous time series data, from lynx trappings
      infstr = 'canlynx.dat' ;
      frmtstr = '%g' ;
      insize = [1,inf] ;
      vicol = [1] ;
      datatran = [] ;
      xstr = 'Year' ;
      ystr = '# Lynx' ;
    elseif iidat == 112 ;
      dtitstr = 'Canadian Lynx - mean' ;
           %  Famous time series data, from lynx trappings
      infstr = 'canlynx.dat' ;
      frmtstr = '%g' ;
      insize = [1,inf] ;
      vicol = [1] ;
      datatran = 'data = data - mean(data) ;' ;
      xstr = 'Year' ;
      ystr = '# Lynx - mean' ;
    elseif iidat == 113 ;
      dtitstr = 'log10(Canadian Lynx + 1)' ;
           %  Famous time series data, from lynx trappings
      infstr = 'canlynx.dat' ;
      frmtstr = '%g' ;
      insize = [1,inf] ;
      vicol = [1] ;
      datatran = 'data = log10(data + 1) ;' ;
      xstr = 'Year' ;
      ystr = 'log10(# Lynx + 1)' ;
    elseif iidat == 114 ;
      dtitstr = 'log10(Canadian Lynx + 1) - mean' ;
           %  Famous time series data, from lynx trappings
      infstr = 'canlynx.dat' ;
      frmtstr = '%g' ;
      insize = [1,inf] ;
      vicol = [1] ;
      datatran = 'data = log10(data + 1) ; data = data - mean(data) ;' ;
      xstr = 'Year' ;
      ystr = 'log10(#L + 1) - mean' ;
    end ;


    disp(['    Loading ' dtitstr ' Data']) ;
      infstr = ['\gauss\steve\data\' infstr] ;
      fid = fopen(infstr,'r') ;
        data = fscanf(fid,frmtstr,insize) ;
        data = data' ;
        %  since data rows are read in as columns
        data = data(:,vicol) ;
      fclose(fid) ;
    disp(['    Read in ' num2str(size(data,1)) ' data points']) ;


    if length(datatran) ~= 0 ;    %  Then do some data transformation
      eval(datatran) ;
    end ;


    gpfoan(data) ;
      subplot(3,1,1) ;
        title(dtitstr) ;
        xlabel(xstr) ;
        ylabel(ystr) ;


    if floor(idat/10) == idat/10 ;    %  then are looping, so pause
      disp(['Any key to continue']) ;
      pause ;
    end ;

  end ;


end ;
