function output = gpvper(data,modlag,lagrange,imode) 
%function [output,moviem,map] = gpvper(data,modlag,lagrange,imode) 
%
% GPVPER, General Purpose Visualization of PERiodicities
%     This does a display for visual insights into periodicities, as:
%         imode = 1:  display seasonal component, deseas'd series,
%                and display of variability explained by other lags
%         imode = 2:  dipsplay "cycle plot", basically a scatterplot
%                over domain [0,lag), with chosen seasonal fit
%         imode = 3;  movie analysis, essentially "morphing through"
%                cycle plot of imode = 2 (recommend choosing a stable r,
%                                         found by imode = 2 analysis)
%
% Inputs:
%     data      -  n x 1 column vector of time series data
%     modlag    -  vector or matrix (optional) of lags whose subspaces
%                     are to be "modded out", before the analysis.
%                     each row (up to 2) corresponds to a subspace.
%                     numbers in a row are lags, pad with zeros
%                     for different numbers in rows.
%                     use a single 0 for no modding out.
%                     WARNING: when this is nonzero, nxn matrices
%                     are created, which can be slow, and 
%                     require a lot of memory.
%     lagrange  -  range of lags to consider (optional)
%                     can be any of:
%                     - single value
%                     - 2 x 1 vector of integer values,
%                            first needs to be <= second,
%                            (only used for imode = 1)
%                     - vector of explicit lags (length ~= 2)
%                            this is most useful for several noninteger
%                            lags.  imode = 2, requires this to have 
%                            length 1,3 or 5
%     imode -      index of mode (optional):
%                     1 - does initial deseasonalization, and display
%                           of variability at a range of lag, together
%                           with deseasonalization at best lag (default)
%                     2 - make "cycle plot", which shows how actual
%                           cycles compare to fit cycles, 
%                           on time scale [0,lag)
%
% Output:
%     graphics only, in a new figure window

%    Copyright (c) J. S. Marron 1998



%  set up parameters
%
n = length(data) ;
maxlag = min(n/2,100) ;
          %  100 is most can comfortaly view, and lags > n/2
          %  are little more interesting than the raw data
if nargin == 1 ; 
          %  only 1 argument input, use defaults

  imodlag = 0 ;
  ilagrange = [2,maxlag] ;
  iimode = 1 ;


  else ;    %  at least 2 arguments input

  imodlag = modlag ;
           %  modlag was specified, use that

  if nargin == 2 ;     %  only 2 arguments input, use defaults for others

    ilagrange = [2,maxlag] ;
    iimode = 1 ;

  else ;    %  at least 3 arguments input

    ilagrange = lagrange ;
           %  lagrange was specified, use that

    if nargin == 3 ;     %  only 3 arguments input, use default iimode
      iimode = 1 ;
    else ;       %  imode was specified, use that
      iimode = imode ;
    end ;

  end ;

end ;




%  calculate preliminary stuff
%
dmin = min(data) ;
dmax = max(data) ;
dmean = mean(data) ;
range = dmax - dmin ;
bottom = dmin - .05 * range ;
top = dmax + .05 * range ;
axhr = .55 * range ;
            %  half range of vertical axis

tss = sum(data.^2) ;
datamm = data - dmean ;
mrss = sum(datamm.^2) ;
df = n - 1 ;
          %  current degrees of freedom in data


if imodlag == 0 ;
  nplotrows = 2 ;
else ;
  nmod = size(imodlag,1) ;
  vmod1 = imodlag(1,:) ;
  vmod1 = vmod1(vmod1 ~= 0) ;
          %  kick out padded zero entries
  nplotrows = 3 ;
  if nmod >= 2 ;    %  then mod out another subspace 
    vmod2 = imodlag(2,:) ;
    vmod2 = vmod2(vmod2 ~= 0) ;
          %  kick out padded zero entries
    nplotrows = 4 ;
  end ;
end ;

if length(ilagrange) == 1 ;
  minlag = ilagrange ;
  maxlag = ilagrange ;
elseif length(ilagrange) == 2 ;
  minlag = ilagrange(1) ;
  maxlag = ilagrange(2) ;
else ;
  minlag = 0 ;
  maxlag = 0 ;
end ;
if maxlag >= n ;
  maxlag = n - 1 ;
end ;
          %  lags n and bigger are not interesting, 
          %  since they return exactly the raw data




fighand = figure ;     %  always create a new figure window and get handle




if iimode == 1 ;    %  then optimize over range of lags

  %  plot raw data 
  %
  subplot(nplotrows,2,1) ;
    plot((1:n)',data,'-',[1; n],[dmean; dmean],'--') ;
      title(['Raw Data (and mean)']) ;
      axis([1, n, bottom, top]) ;
        xt = 1 + .2 * (n - 1) ;
        yt = bottom + .9 * range ;
        tstr = ['Mean is ' num2str(100 * (tss - mrss) / tss) '% of Total SS'] ;
      text(xt, yt, tstr) ;

end ;





if nplotrows > 2 ;  %  then need to subtract off some projections

  %  Get projection matrix onto first subspace
  %
  x = [] ;
  for ilag = 1:length(vmod1) ;  %  loop through lags for first subspace

    %  create columns of x for this lag
    %
    lag = vmod1(ilag) ;
    xl = [] ;
    for iblock = 1:floor(n/lag) ;    %  create size lag identity matrices
      xl = [xl; eye(lag)] ;
    end ;
    nleftover = n - lag * floor(n/lag) ;
    if nleftover > 0 ;
              %  then have some left overs, so need to tack on end part
      xl = [xl; eye(nleftover,lag)] ;
    end ;

    x = [x, xl] ;

  end ;

  ps1 = x * pinv(x' * x) * x' ;
        %  projection matrix for first lagged subspace


  csc = ps1 * datamm ;
        %  project data onto first subspace

  datamm = datamm - csc ;
        %  projection onto orthogonal component

  df = df - trace(ps1) ;
        %  since data now have fewer degrees of freedon

  if iimode == 1 ;    %  then put this into the current plot

    %  Plot Current Seasonal piece
    %
    subplot(nplotrows,2,2) ;
      plot((1:n)',csc) ;
        title(['Seasonal Component, lags ' num2str(vmod1)]) ;
        axis([1, n, -axhr, axhr]) ;
              %  put on same vertical scale for visual comparison
          xt = 1 + .2 * (n - 1) ;
          yt = .8 * axhr ;
          cscss = sum(csc.^2) ;
          tstr = [num2str(100 * cscss / mrss) '% of MR SS'] ;
        text(xt, yt, tstr) ;


    %  Plot residuals from current Seasonal piece
    %
    subplot(nplotrows,2,3) ;
      plot((1:n)',datamm) ;
        title(['Residuals from deseas. lags ' num2str(vmod1)]) ;
        axis([1, n, -axhr, axhr]) ;
              %  put on same vertical scale for visual comparison
          xt = 1 + .2 * (n - 1) ;
          yt = .8 * axhr ;
          dmmss = sum(datamm.^2) ;
          tstr = [num2str(100 * dmmss / mrss) '% of MR SS'] ;
        text(xt, yt, tstr) ;
        hold on ;
          plot([1, n],[0, 0],'--') ;
              %  put dashed line at 0
        hold off ;

  end ;


  if nplotrows == 4 ;    %  then add second decomposition

    %  Get projection matrix onto second subspace
    %
    x = [] ;
    for ilag = 1:length(vmod2) ;  %  loop through lags for first subspace

      %  create columns of x for this lag
      %
      lag = vmod2(ilag) ;
      xl = [] ;
      for iblock = 1:floor(n/lag) ;    %  create size lag idnetity matrices
        xl = [xl; eye(lag)] ;
      end ;
      nleftover = n - lag * floor(n/lag) ;
      if nleftover > 0 ;
              %  then have some left overs, so need to tack on end part
        xl = [xl; eye(nleftover,lag)] ;
      end ;

      x = [x, xl] ;

    end ;


    x = (eye(n) - ps1) * x ;
          %  project design matrix onto orthogonal complement of s1

    ps2ms1 = x * pinv(x' * x) * x' ;
          %  projection matrix for second subspace modulo the first



    csc = ps2ms1 * datamm ;
          %  project data onto second subspace, modulo the first

    datamm = datamm - csc ;
          %  projection onto orthogonal component

    df = df - trace(ps2ms1) ;
          %  since data now have fewer degrees of freedon

    if iimode == 1 ;    %  then add this into the current plot

      %  Plot Current Seasonal piece
      %
      subplot(nplotrows,2,4) ;
        plot((1:n)',csc) ;
          title(['Seasonal Component, lags ' num2str(vmod2)]) ;
          axis([1, n, -axhr, axhr]) ;
                %  put on same vertical scale for visual comparison
          xt = 1 + .2 * (n - 1) ;
          yt = .8 * axhr ;
          cscss = sum(csc.^2) ;
          tstr = [num2str(100 * cscss / mrss) '% of MR SS'] ;
        text(xt, yt, tstr) ;


      %  Plot residuals from current Seasonal piece
      %
      subplot(nplotrows,2,5) ;
        plot((1:n)',datamm) ;
          title(['Residuals from deseas. lags ' num2str(vmod2)]) ;
          axis([1, n, -axhr, axhr]) ;
                %  put on same vertical scale for visual comparison
          xt = 1 + .2 * (n - 1) ;
          yt = .8 * axhr ;
          dmmss = sum(datamm.^2) ;
          tstr = [num2str(100 * dmmss / mrss) '% of MR SS'] ;
        text(xt, yt, tstr) ;
          hold on ;
            plot([1, n],[0, 0],'--') ;
                %  put dashed line at 0
          hold off ;


    end  ;

  end ;

end ;



if iimode == 1 ;    %  then optimize over range of lags

  vperex = [] ;
  bestbetap = 0 ;
  bestdsd = [] ;
  bestsc = [] ;
  bestlag = 0 ;
  if maxlag == 0 ;    %  then use lags in vector ilagrange

    for ilag = 1:length(ilagrange) ;

      lag = ilagrange(ilag) ;

      [dsd,sc] = gpdeseas(datamm,lag) ;
              %  do lag deseasonalization

      perex = 100 * sum((sc - mean(sc)).^2) / mrss ;
              %  percent explained by this lag

      betap = betacdf(perex / 100,lag / 2,(df-lag) / 2) ;
              %  Probability for appropriate beta distribution
     if betap > bestbetap ;     %  then update "best so far" versions
        bestbetap = betap ;
        bestdsd = dsd ;
        bestsc = sc ;
        bestlag = lag ;
      end ;

      vperex = [vperex, perex] ;


    end ;

  else ;    %  then create integer grid of lags, from minlag to maxlag

    for lag = minlag:maxlag ;
  
      [dsd,sc] = gpdeseas(datamm,lag) ;
              %  do lag deseasonalization

      perex = 100 * sum((sc - mean(sc)).^2) / mrss ;
              %  percent explained by this lag

      betap = betacdf(perex / 100,lag / 2,(df-lag) / 2) ;
              %  Probability for appropriate beta distribution

      if betap > bestbetap ;     %  then update "best so far" versions
        bestbetap = betap ;
        bestdsd = dsd ;
        bestsc = sc ;
        bestlag = lag ;
      end ;

      vperex = [vperex, perex] ;

    end ;

  end ;


  %  Plot best Seasonal piece
  %
  subplot(nplotrows,2,2*(nplotrows-2)+2) ;
    plot((1:n)',bestsc) ;
      title(['Best Seasonal Component, lag ' num2str(bestlag)]) ;
      axis([1, n, -axhr, axhr]) ;
            %  put on same vertical scale for visual comparison
        xt = 1 + .2 * (n - 1) ;
        yt = .8 * axhr ;
        bscss = sum(bestsc.^2) ;
        tstr = [num2str(100 * bscss / mrss) '% of MR SS'] ;
      text(xt, yt, tstr) ;


  %  Plot residuals from best Seasonal piece
  %
  subplot(nplotrows,2,2*(nplotrows-2)+3) ;
    plot((1:n)',bestdsd) ;
      title(['Residuals from best deseas. lag ' num2str(bestlag)]) ;
      axis([1, n, -axhr, axhr]) ;
            %  put on same vertical scale for visual comparison
        xt = 1 + .2 * (n - 1) ;
        yt = .8 * axhr ;
        dsdss = sum(bestdsd.^2) ;
        tstr = [num2str(100 * dsdss / mrss) '% of MR SS'] ;
      text(xt, yt, tstr) ;
      hold on ;
        plot([1, n],[0, 0],'--') ;
            %  put dashed line at 0
      hold off ;



  %  Plot Percents Explained
  %
  alpha = 0.05 ;
  if (minlag < maxlag) | (maxlag == 0) ;     %  if there is something to plot
    subplot(nplotrows,2,2*(nplotrows-2)+4) ;
    if maxlag == 0 ;    %  then use lags in vector ilagrange
      plot(ilagrange,vperex) ;

        lowb=betainv(alpha,ilagrange / 2, ...
                               (df-ilagrange) / 2) * 100 ;
        upperb=betainv(1 - alpha,ilagrange / 2, ...
                               (df-ilagrange) / 2) * 100 ;  
            %  pointwise (each end) confidence bands

        nlag = length(ilagrange) ;
        blowb=betainv(alpha/nlag,ilagrange / 2, ...
                               (df-ilagrange) / 2) * 100 ;
        bupperb=betainv(1 - alpha/nlag,ilagrange / 2, ...
                               (df-ilagrange) / 2) * 100 ;      
            %  Bonferroni confidence bands

        hold on ;
          plot(ilagrange,lowb,':');
          plot(ilagrange,upperb,':');

          plot(ilagrange,blowb,'--');
          plot(ilagrange,bupperb,'--');
        hold off ;

        vax = axis ;
        axis([min(ilagrange), max(ilagrange), vax(3), vax(4)]) ;

    else ;    %  then create integer grid of lags, from minlag to maxlag
      plot((minlag:maxlag),vperex) ;

        lowb=betainv(alpha,[minlag:maxlag] / 2, ...
                               (df-[minlag:maxlag]) / 2) * 100 ;
        upperb=betainv(1 - alpha,[minlag:maxlag] / 2, ...
                               (df-[minlag:maxlag]) / 2) * 100 ;  
            %  pointwise (each end) confidence bands

        nlag = maxlag - minlag + 1 ;
        blowb=betainv(alpha/nlag,[minlag:maxlag] / 2, ...
                               (df-[minlag:maxlag]) / 2) * 100 ;
        bupperb=betainv(1 - alpha/nlag,[minlag:maxlag] / 2, ...
                               (df-[minlag:maxlag]) / 2) * 100 ;      
            %  Bonferroni confidence bands

        hold on ;
          plot([minlag:maxlag],lowb,':');
          plot([minlag:maxlag],upperb,':');

          plot([minlag:maxlag],blowb,'--');
          plot([minlag:maxlag],bupperb,'--');
        hold off ;

        vax = axis ;
        axis([minlag, maxlag, vax(3), vax(4)]) ;

    end ;
        title('Percent of SS Explained by these Lags') ;
        hold on ;
          plot([bestlag, bestlag], [vax(3), vax(4)]) ;
            %  overlay vertical line showing best

        hold off ;

  end ;


elseif iimode == 2 ;    %  then make cycle wise plots


  if (length(ilagrange) == 1) | ...
        (length(ilagrange) == 3) | ...
        (length(ilagrange) == 5) ;
            %  then can make cycle plot

    vpmrss = [] ;
          %  vector for keeping percent of MR SS
    for ilag = 1:length(ilagrange) ;

      lag = ilagrange(ilag) ;

      vt = (0:(n-1))' ;
           %  variable representing "time"
      vt = vt / lag ;
           %  maps [0,lag) to [0,1)
      vt = vt - floor(vt) ;
           %  remainder after division
      vt = vt * lag ;
           %  map back to [0,lag)
      bevt = [(vt - lag); vt; (vt + lag)] ;
      bevdata = [datamm; datamm; datamm] ;
           %  boundary extended versions


      [vseasgrid, xgrid] = gpnpr([bevt, bevdata],0,[-lag; (2 * lag)]) ;
           %  seasonal component on boundary extended grid
           %  0 for RSW smooth

      [dsd,sc] = gpdeseas(datamm,lag) ;
              %  do lag deseasonalization, just for Sum of Squares
          bss = sum(sc.^2) ;


      if (length(ilagrange) == 1) ;
        subplot(1,1,1) ;
      elseif (length(ilagrange) == 3) ;
        if ilag == 1 ;
          subplot(2,2,1) ;
        elseif ilag == 2 ;
          subplot(2,2,3) ;
        elseif ilag == 3 ;
          subplot(2,2,4) ;
        end ;
      elseif (length(ilagrange) == 5) ;
        if ilag == 1 ;
          subplot(2,3,1) ;
        elseif ilag == 2 ;
          subplot(2,3,4) ;
        elseif ilag == 3 ;
          subplot(2,3,5) ;
        elseif ilag == 4 ;
          subplot(2,3,6) ;
        elseif ilag == 5 ;
          subplot(2,3,2) ;
        end ;
      end ;

      pmrss = 100 * bss / mrss ;
      vpmrss = [vpmrss; pmrss] ;


      plot(xgrid, vseasgrid) ;
      vachil = get(gca,'Children') ;
        axis([-lag/2,3*lag/2,-axhr,1.5*axhr]) ;
                %  1.5*axhr allows curve showing maxima
        set(vachil(1),'LineWidth',2) ;
        title(['Lag ' num2str(lag) ' Cycle Plot']) ;
          xt = -0.4 * lag ;
          yt = 1.3 * axhr ;
        text(xt,yt,[num2str(pmrss) '% of MR SS']) ;

        hold on ;
          plot([0; 0],[-axhr; 1.5*axhr],'--') ;
          plot([lag; lag],[-axhr; 1.5*axhr],'--') ;

          plvt = vt ;
          plvdata = datamm ;
          plflag = (plvt(2:n) - plvt(1:(n-1)) < 0) ;
               %  has 1s at last points before time comes back
          vmaxloc = [] ;
          while ~isempty(plvt) ;
            [maxplf, ilast] = max(plflag) ;
            if maxplf == 1 ;  %  if there is a next segment
              tpvt = plvt(1:ilast) ;
              tpvdata = plvdata(1:ilast) ;
            else ;      %  then doing last piece
              tpvt = plvt ;
              tpvdata = plvdata ;
            end ;

            tpn = length(tpvt) ;
            if tpn > 1 ;   %  need at least 2 points to talk about max
                           %  since can get third from wrap
              [temp,imax] = max(tpvdata) ;
              if imax == 1 ;    %  then max on leftend, so wrap
                qcoeffs = polyfit([tpvt(tpn)-lag;tpvt(1:2)],...
                                       [tpvdata(tpn);tpvdata(1:2)],2) ;
              elseif imax == tpn ;    %  then max on rightend, so wrap
                qcoeffs = polyfit([tpvt(tpn-1:tpn);tpvt(1)+lag],...
                                       [tpvdata(tpn-1:tpn);tpvdata(1)],2) ;
              else ;    %  then max in interior
                qcoeffs = polyfit(tpvt((imax-1):(imax+1)),...
                                       tpvdata((imax-1):(imax+1)),2) ;
              end ;
              vmaxloc = [vmaxloc; -qcoeffs(2) / (2* qcoeffs(1))] ;
                  %  location of max, for overlaying later
                  %      (uses -b/2a formula)
            end ;

            tpvt = [(tpvt-lag); tpvt; (tpvt+lag)] ;
            tpvdata = [tpvdata; tpvdata; tpvdata] ;
            plot(tpvt,tpvdata,':') ;

            if maxplf == 1 ;  %  if there is a next segment
              plvt = plvt((ilast+1):length(plvt)) ;
              plvdata = plvdata((ilast+1):length(plvdata)) ;
              plflag = plflag((ilast+1):length(plflag)) ;
            else ;      %  then doing last piece
              plvt = [] ;
              plvdata = [] ;
              plflag = [] ;
            end ;
          end ;

          %  add plot of max locations
          %
          nmaxt = length(vmaxloc) ;
          for im = 2:nmaxt ;     %  loop through maxes 
            if vmaxloc(im) < vmaxloc(im - 1) - lag/2 ; 
                    %  then goes too far down, 
                    %  so move rest of sequence to right
              vmaxloc(im:nmaxt) = vmaxloc(im:nmaxt) + lag ;
            elseif vmaxloc(im) > vmaxloc(im - 1) + lag/2 ; 
                    %  then goes too far up, 
                    %  so move rest of sequence to left
              vmaxloc(im:nmaxt) = vmaxloc(im:nmaxt) - lag ;
            end ;
          end ;
          vmaxt = (1:nmaxt)' ;
          vmaxt = (1 + 0.5 * (vmaxt / nmaxt)) * axhr ;
          plot(vmaxloc - lag,vmaxt,'-c') ;
          plot(vmaxloc,vmaxt,'-c') ;
          plot(vmaxloc + lag,vmaxt,'-c') ;

        hold off ;

    end ;


    if (length(ilagrange) == 3) ;
      subplot(2,2,2) ;
        bar((1:3)',vpmrss) ;
        axis([0.3,3.7,0,102]) ;
        title('% of MR SS, for each pic') ;

    elseif (length(ilagrange) == 5) ;
  
    end ;


  else ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from gpvert.m:       !!!') ;
    disp('!!!       Number of input lags   !!!') ;
    disp('!!!       must be 1, 3 or 5      !!!') ;
    disp('!!!       No plot made           !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

  end ;


elseif iimode == 3 ;   %  then make movies, of cycle wise plots

  if (length(ilagrange) == 1) ;   %  then can make movie of cycle plot

    disp('constructing movie') ;
    
    
    set(fighand,'position',[240,300,280,210]) ;
 
    
    vpmrss = [] ;
          %  vector for keeping percent of MR SS

    lag = ilagrange(1) ;


    %  first get overall seasonal part to overalay
    %
    vt = (0:(n-1))' ;
          %  variable representing "time"
    vt = vt / lag ;
          %  maps [0,lag) to [0,1)
    vt = vt - floor(vt) ;
          %  remainder after division
    vt = vt * lag ;
          %  map back to [0,lag)
    bevt = [(vt - lag); vt; (vt + lag)] ;
    bevdata = [datamm; datamm; datamm] ;
          %  boundary extended versions

    [vseasgrid, xgrid] = gpnpr([bevt, bevdata],0,[-lag; (2 * lag)]) ;
          %  seasonal component on boundary extended grid
          %  0 for RSW smooth


    %  get flag for cutting time series into cycles
    %
    plvt = vt ;
    plvdata = datamm ;
    plflag = (plvt(2:n) - plvt(1:(n-1)) < 0) ;
          %  has 1s at last points before time comes back

    %  cut off stuff after last 1
    %
    [temp,indlast1] = max(flipud(plflag)) ;
    indlast1 = length(plflag) + 1 - indlast1 ;
          %  index of last 1 in plflag
    plvt = plvt(1:indlast1) ;
    plvdata = plvdata(1:indlast1) ;
    plflag = plflag(1:indlast1) ;
    

    nmaxt = sum(plflag) ;
    if (nmaxt > 5) ;    %  Then have at least 5 cycles, so make movie


      % first put down plot basics
      %
      subplot(1,1,1) ;
      plot(xgrid, vseasgrid,'m-') ;
      vachil = get(gca,'Children') ;
        axis([-lag/2,3*lag/2,-axhr,1.5*axhr]) ;
                %  1.5*axhr allows curve showing maxima
        set(vachil(1),'LineWidth',4) ;
        title(['Lag ' num2str(lag) ' Cycle Plot']) ;

        hold on ;
          plot([0; 0],[-axhr; 1.5*axhr],'b--') ;
          plot([lag; lag],[-axhr; 1.5*axhr],'b--') ;
        hold off ;

     


      %  loop through first 5 cycles
      %
      vmaxloc = [] ;
      th=[];
      tmaxloc=[]; 
      for icyc = 1:5 ;

        [maxplf, ilast] = max(plflag) ;
            %  ilast is indx of first one (i.e. location of end of cycle)
        tpvt = plvt(1:ilast) ;
            %  current vector of times
        tpvdata = plvdata(1:ilast) ;
            %  current vector of heights
        plvt = plvt((ilast+1):length(plvt)) ;
            %  remaining times
        plvdata = plvdata((ilast+1):length(plvdata)) ;
            %  remaining heights
        plflag = plflag((ilast+1):length(plflag)) ;
            %  remaining flags


        %  get location of current maximum
        %
        tpn = length(tpvt) ;
        if tpn > 1 ;   %  need at least 2 points to talk about max
                       %  since can get third from wrap
          [temp,imax] = max(tpvdata) ;

          txcs = [cos(2*pi*[1:lag]'/lag),sin(2*pi*[1:lag]'/lag)];
            % design matrix: txcs=[cos(omega*t),sin(omega*t)]

          tab =inv(txcs'*txcs)*txcs'*tpvdata;
            % estimated coefficients of a and b 

          th = [th;sqrt(sum(tab.^2))];  
            % height of this cycle presented in Fourier frequency;

          tmaxloc = [tmaxloc;(pi/2-atan(tab(1)/tab(2)))*lag/(2*pi)];
            % maximum location presented in Fourier frequency;
          
          if imax == 1 ;    %  then max on leftend, so wrap
            qcoeffs = polyfit([tpvt(tpn)-lag;tpvt(1:2)],...
                                   [tpvdata(tpn);tpvdata(1:2)],2) ;
          elseif imax == tpn ;    %  then max on rightend, so wrap
            qcoeffs = polyfit([tpvt(tpn-1:tpn);tpvt(1)+lag],...
                                   [tpvdata(tpn-1:tpn);tpvdata(1)],2) ;
          else ;    %  then max in interior
            qcoeffs = polyfit(tpvt((imax-1):(imax+1)),...
                                   tpvdata((imax-1):(imax+1)),2) ;
          end ;
          vmaxloc = [vmaxloc; -qcoeffs(2) / (2* qcoeffs(1))] ;
              %  location of max, for overlaying later
              %      (uses -b/2a formula)
        end ;


        %  replicate the data and plot
        %
        tpvt = [(tpvt-lag); tpvt; (tpvt+lag)] ;
        tpvdata = [tpvdata; tpvdata; tpvdata] ;
        hold on ;
          plot(tpvt,tpvdata,'g-') ;
          vachil = get(gca,'Children') ;
          if icyc == 2 | icyc == 4 ;
            set(vachil(1),'LineWidth',2) ;
          elseif icyc == 3 ;
            set(vachil(1),'LineWidth',4) ;
          end ;
        hold off ;


        eval(['tpvt' num2str(icyc) ' = tpvt ;']) ;
        eval(['tpvdata' num2str(icyc) ' = tpvdata ;']) ;

      end ;


      %  add plot of max locations
      %
      hold on ;
        for im = 2:5 ;     %  loop through maxes 
          if vmaxloc(im) < vmaxloc(im - 1) - lag/2 ; 
                  %  then goes too far down, 
                  %  so move rest of sequence to right
            vmaxloc(im:length(vmaxloc)) = vmaxloc(im:length(vmaxloc)) + lag ;
          elseif vmaxloc(im) > vmaxloc(im - 1) + lag/2 ; 
                  %  then goes too far up, 
                  %  so move rest of sequence to left
            vmaxloc(im:length(vmaxloc)) = vmaxloc(im:length(vmaxloc)) - lag ;
          end ;
        end ;
%        vmaxt = (1:length(vmaxloc))' ;
%        vmaxt = (1 + 0.5 * (vmaxt / nmaxt)) * axhr ;
%        plot(vmaxloc - lag,vmaxt,'-c') ;
%        plot(vmaxloc(2) - lag,vmaxt(2),'co') ;
%        plot(vmaxloc(3) - lag,vmaxt(3),'co') ;
%        plot(vmaxloc(4) - lag,vmaxt(4),'co') ;
%        vachil = get(gca,'Children') ;
%        set(vachil(1),'MarkerSize',3) ;
%        set(vachil(3),'MarkerSize',3) ;


%        plot(vmaxloc,vmaxt,'-c') ;
%        plot(vmaxloc(2),vmaxt(2),'co') ;
%        plot(vmaxloc(3),vmaxt(3),'co') ;
%        plot(vmaxloc(4),vmaxt(4),'co') ;
%        vachil = get(gca,'Children') ;
%        set(vachil(1),'MarkerSize',2) ;
%        set(vachil(3),'MarkerSize',2) ;

%        plot(vmaxloc + lag,vmaxt,'-c') ;
%        plot(vmaxloc(2) + lag,vmaxt(2),'co') ;
%        plot(vmaxloc(3) + lag,vmaxt(3),'co') ;
%        plot(vmaxloc(4) + lag,vmaxt(4),'co') ;
%        vachil = get(gca,'Children') ;
%        set(vachil(1),'MarkerSize',2) ;
%        set(vachil(3),'MarkerSize',2) ;
        

        vmaxt = (1:length(tmaxloc))' ;
        vmaxt = (1 + 0.5 * (vmaxt / nmaxt)) * axhr ;
        plot(tmaxloc - lag,vmaxt,'-c') ;
        plot(tmaxloc(2) - lag,vmaxt(2),'co') ;
        plot(tmaxloc(3) - lag,vmaxt(3),'co') ;
        plot(tmaxloc(4) - lag,vmaxt(4),'co') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',3) ;
        set(vachil(3),'MarkerSize',3) ;


        plot(tmaxloc,vmaxt,'-c') ;
        plot(tmaxloc(2),vmaxt(2),'co') ;
        plot(tmaxloc(3),vmaxt(3),'co') ;
        plot(tmaxloc(4),vmaxt(4),'co') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',2) ;
        set(vachil(3),'MarkerSize',2) ;

        plot(tmaxloc + lag,vmaxt,'-c') ;
        plot(tmaxloc(2) + lag,vmaxt(2),'co') ;
        plot(tmaxloc(3) + lag,vmaxt(3),'co') ;
        plot(tmaxloc(4) + lag,vmaxt(4),'co') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',2) ;
        set(vachil(3),'MarkerSize',2) ;        
      hold off ;



      moviem = [getframe(gcf) getframe(gcf) getframe(gcf) getframe(gcf) getframe(gcf) ];
          %  first frame of movie



 
      %  loop through next frames of movie
      %
      while ~isempty(plvt) ;
        [maxplf, ilast] = max(plflag) ;
        if maxplf == 1 ;  %  if there is a next segment
          tpvt = plvt(1:ilast) ;
          tpvdata = plvdata(1:ilast) ;





        else ;      %  then doing last piece
          tpvt = plvt ;
          tpvdata = plvdata ;
        end ;

        %  get location of current maximum
        %
        tpn = length(tpvt) ;
        if tpn > 1 ;   %  need at least 2 points to talk about max
                       %  since can get third from wrap
          [temp,imax] = max(tpvdata) ;
         txcs = [cos(2*pi*[1:lag]'/lag),sin(2*pi*[1:lag]'/lag)];
            % design matrix: txcs=[cos(omega*t),sin(omega*t)]

          tab =inv(txcs'*txcs)*txcs'*tpvdata;
            % estimated coefficients of a and b 

          th = [th;sqrt(sum(tab.^2))];  
            % height of this cycle presented in Fourier frequency;

          tmaxloc = [tmaxloc;(pi/2-atan(tab(1)/tab(2)))*lag/(2*pi)];
            % maximum location presented in Fourier frequency;
             
          if imax == 1 ;    %  then max on leftend, so wrap
            qcoeffs = polyfit([tpvt(tpn)-lag;tpvt(1:2)],...
                                   [tpvdata(tpn);tpvdata(1:2)],2) ;
          elseif imax == tpn ;    %  then max on rightend, so wrap
            qcoeffs = polyfit([tpvt(tpn-1:tpn);tpvt(1)+lag],...
                                   [tpvdata(tpn-1:tpn);tpvdata(1)],2) ;
          else ;    %  then max in interior
            qcoeffs = polyfit(tpvt((imax-1):(imax+1)),...
                                   tpvdata((imax-1):(imax+1)),2) ;
          end ;
          vmaxloc = [vmaxloc; -qcoeffs(2) / (2* qcoeffs(1))] ;
              %  location of max, for overlaying later
              %      (uses -b/2a formula)
        end ;


        %  adjust max location
        if vmaxloc(length(vmaxloc)) < ...
                                 vmaxloc(length(vmaxloc) - 1) - lag/2 ; 
                %  then goes too far down, 
                %  so move to right
          vmaxloc(length(vmaxloc)) = vmaxloc(length(vmaxloc)) + lag ;
        elseif vmaxloc(length(vmaxloc)) > ...
                                 vmaxloc(length(vmaxloc) - 1) + lag/2 ; 
                %  then goes too far up, 
                %  so move to left
          vmaxloc(length(vmaxloc)) = vmaxloc(length(vmaxloc)) - lag ;
        end ;

        tpvt = [(tpvt-lag); tpvt; (tpvt+lag)] ;
        tpvdata = [tpvdata; tpvdata; tpvdata] ;

        tpvt1 = tpvt2 ;
        tpvt2 = tpvt3 ;
        tpvt3 = tpvt4 ;
        tpvt4 = tpvt5 ;
        tpvt5 = tpvt ;

        tpvdata1 = tpvdata2 ;
        tpvdata2 = tpvdata3 ;
        tpvdata3 = tpvdata4 ;
        tpvdata4 = tpvdata5 ;
        tpvdata5 = tpvdata ;

        %  make next movie frame
        %
        subplot(1,1,1) ;
        plot(xgrid, vseasgrid,'m-') ;
        vachil = get(gca,'Children') ;
          axis([-lag/2,3*lag/2,-axhr,1.5*axhr]) ;
                  %  1.5*axhr allows curve showing maxima
          set(vachil(1),'LineWidth',4) ;
          title(['Lag ' num2str(lag) ' Cycle Plot']) ;
  
          hold on ;
            plot([0; 0],[-axhr; 1.5*axhr],'b--') ;
            plot([lag; lag],[-axhr; 1.5*axhr],'b--') ;

            plot(tpvt1,tpvdata1,'g-') ;
            plot(tpvt2,tpvdata2,'g-') ;
            plot(tpvt3,tpvdata3,'g-') ;
            plot(tpvt4,tpvdata4,'g-') ;
            plot(tpvt5,tpvdata5,'g-') ;
            vachil = get(gca,'Children') ;
            set(vachil(2),'LineWidth',2) ;
            set(vachil(3),'LineWidth',4) ;
            set(vachil(4),'LineWidth',2) ;

  
            %  add plot of max locations
            %
%            vmaxt = (1:length(vmaxloc))' ;
%            vmaxt = (1 + 0.5 * (vmaxt / nmaxt)) * axhr ;
%            plot(vmaxloc - lag,vmaxt,'-c') ;
%            plot(vmaxloc(length(vmaxt)-1) - lag,vmaxt(length(vmaxt)-1),'co') ;
%            plot(vmaxloc(length(vmaxt)-2) - lag,vmaxt(length(vmaxt)-2),'co') ;
%            plot(vmaxloc(length(vmaxt)-3) - lag,vmaxt(length(vmaxt)-3),'co') ;
%            vachil = get(gca,'Children') ;
%            set(vachil(1),'MarkerSize',3) ;
%            set(vachil(3),'MarkerSize',3) ;

%            plot(vmaxloc,vmaxt,'-c') ;
%            plot(vmaxloc(length(vmaxt)-1),vmaxt(length(vmaxt)-1),'co') ;
%            plot(vmaxloc(length(vmaxt)-2),vmaxt(length(vmaxt)-2),'co') ;
%            plot(vmaxloc(length(vmaxt)-3),vmaxt(length(vmaxt)-3),'co') ;
%            vachil = get(gca,'Children') ;
%            set(vachil(1),'MarkerSize',3) ;
%            set(vachil(3),'MarkerSize',3) ;

%            plot(vmaxloc + lag,vmaxt,'-c') ;
%            plot(vmaxloc(length(vmaxt)-1) + lag,vmaxt(length(vmaxt)-1),'co') ;
%            plot(vmaxloc(length(vmaxt)-2) + lag,vmaxt(length(vmaxt)-2),'co') ;
%            plot(vmaxloc(length(vmaxt)-3) + lag,vmaxt(length(vmaxt)-3),'co') ;
%            vachil = get(gca,'Children') ;
%            set(vachil(1),'MarkerSize',3) ;
%            set(vachil(3),'MarkerSize',3) ;
%            
            vmaxt = (1:length(tmaxloc))' ;
            vmaxt = (1 + 0.5 * (vmaxt / nmaxt)) * axhr ;
            plot(tmaxloc - lag,vmaxt,'-c') ;
            plot(tmaxloc(length(vmaxt)-1) - lag,vmaxt(length(vmaxt)-1),'co') ;
            plot(tmaxloc(length(vmaxt)-2) - lag,vmaxt(length(vmaxt)-2),'co') ;
            plot(tmaxloc(length(vmaxt)-3) - lag,vmaxt(length(vmaxt)-3),'co') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;

            plot(tmaxloc,vmaxt,'-c') ;
            plot(tmaxloc(length(vmaxt)-1),vmaxt(length(vmaxt)-1),'co') ;
            plot(tmaxloc(length(vmaxt)-2),vmaxt(length(vmaxt)-2),'co') ;
            plot(tmaxloc(length(vmaxt)-3),vmaxt(length(vmaxt)-3),'co') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;

            plot(tmaxloc + lag,vmaxt,'-c') ;
            plot(tmaxloc(length(vmaxt)-1) + lag,vmaxt(length(vmaxt)-1),'co') ;
            plot(tmaxloc(length(vmaxt)-2) + lag,vmaxt(length(vmaxt)-2),'co') ;
            plot(tmaxloc(length(vmaxt)-3) + lag,vmaxt(length(vmaxt)-3),'co') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;
            
  
          hold off ;





        moviem = [moviem,getframe(gcf),getframe(gcf),...
                     getframe(gcf),getframe(gcf),getframe(gcf)] ;



        if maxplf == 1 ;  %  if there is a next segment
          plvt = plvt((ilast+1):length(plvt)) ;
          plvdata = plvdata((ilast+1):length(plvdata)) ;
          plflag = plflag((ilast+1):length(plflag)) ;
        else ;      %  then doing last piece
          plvt = [] ;
          plvdata = [] ;
          plflag = [] ;
        end ;
  

      end ;


      
      moviem=[moviem,moviem(:,[size(moviem,2):(-1):1])];
      
      movie(fighand,moviem,-2) ;
            %  plays movie 2 times, back and forth
 
      map=colormap;
      mpgwrite(moviem,map,'D411c1new.mpg')
      output=[tmaxloc,th];

    else ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Error from gpvert.m:           !!!') ;
      disp('!!!       Number of cycles < 5       !!!') ;
      disp('!!!       No movie made              !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

    end ;


   end;

 

elseif iimode == 4 ;   %  then  make scatter plot of height vs max-location

  if (length(ilagrange) == 1) ;   %  then make scatter plot of height vs max-loc

   set(fighand,'position',[240,300,280,210]) ;
   
    vpmrss = [] ;
          %  vector for keeping percent of MR SS

    lag = ilagrange(1) ;


    %  first get overall seasonal part to overalay
    %
    vt = (0:(n-1))' ;
          %  variable representing "time"
    vt = vt / lag ;
          %  maps [0,lag) to [0,1)
    vt = vt - floor(vt) ;
          %  remainder after division
    vt = vt * lag ;
          %  map back to [0,lag)
    bevt = [(vt - lag); vt; (vt + lag)] ;
    bevdata = [datamm; datamm; datamm] ;
          %  boundary extended versions

    [vseasgrid, xgrid] = gpnpr([bevt, bevdata],0,[-lag; (2 * lag)]) ;
          %  seasonal component on boundary extended grid
          %  0 for RSW smooth


    %  get flag for cutting time series into cycles
    %
    plvt = vt ;
    plvdata = datamm ;
    plflag = (plvt(2:n) - plvt(1:(n-1)) < 0) ;
          %  has 1s at last points before time comes back

    %  cut off stuff after last 1
    %
    [temp,indlast1] = max(flipud(plflag)) ;
    indlast1 = length(plflag) + 1 - indlast1 ;
          %  index of last 1 in plflag
    plvt = plvt(1:indlast1) ;
    plvdata = plvdata(1:indlast1) ;
    plflag = plflag(1:indlast1) ;
    

    nmaxt = sum(plflag) ;
    if (nmaxt > 5) ;    %  Then have at least 5 cycles, so make movie


      % first put down plot basics
      %
      subplot(1,1,1) ;
      plot(xgrid, vseasgrid,'m-') ;
      vachil = get(gca,'Children') ;
        axis([-lag/2,3*lag/2,-axhr,1.5*axhr]) ;
                %  1.5*axhr allows curve showing maxima
        set(vachil(1),'LineWidth',4) ;
        title(['Lag ' num2str(lag) ' Cycle Plot']) ;

        hold on ;
          plot([0; 0],[-axhr; 1.5*axhr],'b--') ;
          plot([lag; lag],[-axhr; 1.5*axhr],'b--') ;
        hold off ;

     


      %  loop through first 5 cycles
      %
      vmaxloc = [] ;
      height = [];
      frameno=0;
      iframe=[];
      for icyc = 1:5 ;

        [maxplf, ilast] = max(plflag) ;
            %  ilast is indx of first one (i.e. location of end of cycle)
        tpvt = plvt(1:ilast) ;
            %  current vector of times
        tpvdata = plvdata(1:ilast) ;
            %  current vector of heights
        plvt = plvt((ilast+1):length(plvt)) ;
            %  remaining times
        plvdata = plvdata((ilast+1):length(plvdata)) ;
            %  remaining heights
        plflag = plflag((ilast+1):length(plflag)) ;
            %  remaining flags


        %  get location of current maximum
        %
        tpn = length(tpvt) ;
        if tpn > 1 ;   %  need at least 2 points to talk about max
                       %  since can get third from wrap
          [temp,imax] = max(tpvdata) ;
          if imax == 1 ;    %  then max on leftend, so wrap
            qcoeffs = polyfit([tpvt(tpn)-lag;tpvt(1:2)],...
                                   [tpvdata(tpn);tpvdata(1:2)],2) ;
          elseif imax == tpn ;    %  then max on rightend, so wrap
            qcoeffs = polyfit([tpvt(tpn-1:tpn);tpvt(1)+lag],...
                                   [tpvdata(tpn-1:tpn);tpvdata(1)],2) ;
          else ;    %  then max in interior
            qcoeffs = polyfit(tpvt((imax-1):(imax+1)),...
                                   tpvdata((imax-1):(imax+1)),2) ;
          end ;
          vmaxloc = [vmaxloc; -qcoeffs(2) / (2* qcoeffs(1))] ;
              %  location of max, for overlaying later
              %      (uses -b/2a formula)
         height = [height; qcoeffs(3) - qcoeffs(2)^2 /(4 * qcoeffs(1))];
              % height corresponding the max-location
              %   (uses c - b^2 /(4a) )  
         frameno = frameno + 1;    % frame number
         iframe = [iframe; frameno];                 
        end ;


        %  replicate the data and plot
        %
        tpvt = [(tpvt-lag); tpvt; (tpvt+lag)] ;
        tpvdata = [tpvdata; tpvdata; tpvdata] ;
        hold on ;
          plot(tpvt,tpvdata,'g-') ;
          vachil = get(gca,'Children') ;
          if icyc == 2 | icyc == 4 ;
            set(vachil(1),'LineWidth',2) ;
          elseif icyc == 3 ;
            set(vachil(1),'LineWidth',4) ;
          end ;
        hold off ;


        eval(['tpvt' num2str(icyc) ' = tpvt ;']) ;
        eval(['tpvdata' num2str(icyc) ' = tpvdata ;']) ;

      end ;


      %  add plot of max locations
      %
      hold on ;
        for im = 2:5 ;     %  loop through maxes 
          if vmaxloc(im) < vmaxloc(im - 1) - lag/2 ; 
                  %  then goes too far down, 
                  %  so move rest of sequence to right
            vmaxloc(im:length(vmaxloc)) = vmaxloc(im:length(vmaxloc)) + lag ;
          elseif vmaxloc(im) > vmaxloc(im - 1) + lag/2 ; 
                  %  then goes too far up, 
                  %  so move rest of sequence to left
            vmaxloc(im:length(vmaxloc)) = vmaxloc(im:length(vmaxloc)) - lag ;
          end ;
        end ;
        vmaxt = (1:length(vmaxloc))' ;
        vmaxt = (1 + 0.5 * (vmaxt / nmaxt)) * axhr ;
        plot(vmaxloc - lag,vmaxt,'-c') ;
        plot(vmaxloc(2) - lag,vmaxt(2),'co') ;
        plot(vmaxloc(3) - lag,vmaxt(3),'co') ;
        plot(vmaxloc(4) - lag,vmaxt(4),'co') ;
        vachil = get(gca,'Children') ; 
        set(vachil(1),'MarkerSize',3) ;
        set(vachil(3),'MarkerSize',3) ;

        plot(vmaxloc - lag,height,':r') ;
        plot(vmaxloc(2) - lag,height(2),'ro') ;
        plot(vmaxloc(3) - lag,height(3),'r<') ;
        plot(vmaxloc(4) - lag,height(4),'ro') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',2) ;
        set(vachil(3),'MarkerSize',2) ;        
      
        plot(vmaxloc,vmaxt,'-c') ;
        plot(vmaxloc(2),vmaxt(2),'co') ;
        plot(vmaxloc(3),vmaxt(3),'co') ;
        plot(vmaxloc(4),vmaxt(4),'co') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',2) ;
        set(vachil(3),'MarkerSize',2) ;
  
        plot(vmaxloc, height,':b') ;
        plot(vmaxloc(2), height(2),'bo') ;
        plot(vmaxloc(3),height(3),'b.') ;
        plot(vmaxloc(4),height(4),'bo') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',3) ;
        set(vachil(3),'MarkerSize',3) ;  
        
        plot(vmaxloc + lag,vmaxt,'-c') ;
        plot(vmaxloc(2) + lag,vmaxt(2),'co') ;
        plot(vmaxloc(3) + lag,vmaxt(3),'co') ;
        plot(vmaxloc(4) + lag,vmaxt(4),'co') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',2) ;
        set(vachil(3),'MarkerSize',2) ;

        plot(vmaxloc + lag,height,':r') ;
        plot(vmaxloc(2) + lag,height(2),'ro') ;
        plot(vmaxloc(3) + lag,height(3),'r.') ;
        plot(vmaxloc(4) + lag,height(4),'ro') ;
        vachil = get(gca,'Children') ;
        set(vachil(1),'MarkerSize',3) ;
        set(vachil(3),'MarkerSize',3) ;    
      
      hold off ;



      moviem = [getframe(gcf) getframe(gcf) getframe(gcf) getframe(gcf) getframe(gcf) ];
          %  first frame of movie



 
      %  loop through next frames of movie
      %
      while ~isempty(plvt) ;
        [maxplf, ilast] = max(plflag) ;
        if maxplf == 1 ;  %  if there is a next segment
          tpvt = plvt(1:ilast) ;
          tpvdata = plvdata(1:ilast) ;





        else ;      %  then doing last piece
          tpvt = plvt ;
          tpvdata = plvdata ;
        end ;

        %  get location of current maximum
        %
        tpn = length(tpvt) ;
        if tpn > 1 ;   %  need at least 2 points to talk about max
                       %  since can get third from wrap
          [temp,imax] = max(tpvdata) ;
          if imax == 1 ;    %  then max on leftend, so wrap
            qcoeffs = polyfit([tpvt(tpn)-lag;tpvt(1:2)],...
                                   [tpvdata(tpn);tpvdata(1:2)],2) ;
          elseif imax == tpn ;    %  then max on rightend, so wrap
            qcoeffs = polyfit([tpvt(tpn-1:tpn);tpvt(1)+lag],...
                                   [tpvdata(tpn-1:tpn);tpvdata(1)],2) ;
          else ;    %  then max in interior
            qcoeffs = polyfit(tpvt((imax-1):(imax+1)),...
                                   tpvdata((imax-1):(imax+1)),2) ;
          end ;
          vmaxloc = [vmaxloc; -qcoeffs(2) / (2* qcoeffs(1))] ;
              %  location of max, for overlaying later
              %      (uses -b/2a formula)
          height = [height; qcoeffs(3) - qcoeffs(2)^2 /(4 * qcoeffs(1))];
              % height corresponding the max-location
              %   (uses c - b^2 /(4a) )
         frameno = frameno + 1;    % frame number
         iframe = [iframe; frameno];     
        end ;


        %  adjust max location
        if vmaxloc(length(vmaxloc)) < ...
                                 vmaxloc(length(vmaxloc) - 1) - lag/2 ; 
                %  then goes too far down, 
                %  so move to right
          vmaxloc(length(vmaxloc)) = vmaxloc(length(vmaxloc)) + lag ;
        elseif vmaxloc(length(vmaxloc)) > ...
                                 vmaxloc(length(vmaxloc) - 1) + lag/2 ; 
                %  then goes too far up, 
                %  so move to left
          vmaxloc(length(vmaxloc)) = vmaxloc(length(vmaxloc)) - lag ;
        end ;

        tpvt = [(tpvt-lag); tpvt; (tpvt+lag)] ;
        tpvdata = [tpvdata; tpvdata; tpvdata] ;

        tpvt1 = tpvt2 ;
        tpvt2 = tpvt3 ;
        tpvt3 = tpvt4 ;
        tpvt4 = tpvt5 ;
        tpvt5 = tpvt ;

        tpvdata1 = tpvdata2 ;
        tpvdata2 = tpvdata3 ;
        tpvdata3 = tpvdata4 ;
        tpvdata4 = tpvdata5 ;
        tpvdata5 = tpvdata ;

        %  make next movie frame
        %
        subplot(1,1,1) ;
        plot(xgrid, vseasgrid,'m-') ;
        vachil = get(gca,'Children') ;
          axis([-lag/2,3*lag/2,-axhr,1.5*axhr]) ;
                  %  1.5*axhr allows curve showing maxima
          set(vachil(1),'LineWidth',4) ;
          title(['Lag ' num2str(lag) ' Cycle Plot']) ;
  
          hold on ;
            plot([0; 0],[-axhr; 1.5*axhr],'b--') ;
            plot([lag; lag],[-axhr; 1.5*axhr],'b--') ;

            plot(tpvt1,tpvdata1,'g-') ;
            plot(tpvt2,tpvdata2,'g-') ;
            plot(tpvt3,tpvdata3,'g-') ;
            plot(tpvt4,tpvdata4,'g-') ;
            plot(tpvt5,tpvdata5,'g-') ;
            vachil = get(gca,'Children') ;
            set(vachil(2),'LineWidth',2) ;
            set(vachil(3),'LineWidth',4) ;
            set(vachil(4),'LineWidth',2) ;

  
            %  add plot of max locations
            %
            vmaxt = (1:length(vmaxloc))' ;
            vmaxt = (1 + 0.5 * (vmaxt / nmaxt)) * axhr ;
            plot(vmaxloc - lag,vmaxt,'-c') ;
            plot(vmaxloc(length(vmaxt)-1) - lag,vmaxt(length(vmaxt)-1),'co') ;
            plot(vmaxloc(length(vmaxt)-2) - lag,vmaxt(length(vmaxt)-2),'co') ;
            plot(vmaxloc(length(vmaxt)-3) - lag,vmaxt(length(vmaxt)-3),'co') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;
            
            plot(vmaxloc - lag,height,':r') ;
            plot(vmaxloc(length(vmaxt)-1) - lag,height(length(vmaxt)-1),'ro') ;
            plot(vmaxloc(length(vmaxt)-2) - lag,height(length(vmaxt)-1),'r<') ;
            plot(vmaxloc(length(vmaxt)-1) - lag,height(length(vmaxt)-1),'ro') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;             

            plot(vmaxloc,vmaxt,'-c') ;
            plot(vmaxloc(length(vmaxt)-1),vmaxt(length(vmaxt)-1),'co') ;
            plot(vmaxloc(length(vmaxt)-2),vmaxt(length(vmaxt)-2),'co') ;
            plot(vmaxloc(length(vmaxt)-3),vmaxt(length(vmaxt)-3),'co') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;
            
            plot(vmaxloc,height,':b') ;
            plot(vmaxloc(length(vmaxt)-1),height(length(vmaxt)-1),'bo') ;
            plot(vmaxloc(length(vmaxt)-2),height(length(vmaxt)-2),'b<') ;
            plot(vmaxloc(length(vmaxt)-3),height(length(vmaxt)-3),'bo') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;   
            
            plot(vmaxloc + lag,vmaxt,'-c') ;
            plot(vmaxloc(length(vmaxt)-1) + lag,vmaxt(length(vmaxt)-1),'co') ;
            plot(vmaxloc(length(vmaxt)-2) + lag,vmaxt(length(vmaxt)-2),'co') ;
            plot(vmaxloc(length(vmaxt)-3) + lag,vmaxt(length(vmaxt)-3),'co') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;
            
            plot(vmaxloc + lag,height,':r') ;
            plot(vmaxloc(length(vmaxt)-1) + lag,height(length(vmaxt)-1),'ro') ;
            plot(vmaxloc(length(vmaxt)-2) + lag,height(length(vmaxt)-2),'r<') ;
            plot(vmaxloc(length(vmaxt)-3) + lag,height(length(vmaxt)-3),'ro') ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'MarkerSize',3) ;
            set(vachil(3),'MarkerSize',3) ;            
  
          hold off ;




      moviem = [moviem,getframe(gcf),getframe(gcf),...
                     getframe(gcf),getframe(gcf),getframe(gcf)] ;



        if maxplf == 1 ;  %  if there is a next segment
          plvt = plvt((ilast+1):length(plvt)) ;
          plvdata = plvdata((ilast+1):length(plvdata)) ;
          plflag = plflag((ilast+1):length(plflag)) ;
        else ;      %  then doing last piece
          plvt = [] ;
          plvdata = [] ;
          plflag = [] ;
        end ;
  

      end ;

 
   

      moviem=[moviem,moviem(:,[size(moviem,2):(-1):1])];
      
      movie(fighand,moviem,-2) ;
            %  plays movie 2 times, back and forth
      output=[iframe, vmaxloc];
%      output=[iframe, vmaxloc, height,vmaxt];
      map=colormap;
%      mpgwrite(moviem,map,'test.mpg')
      

    else ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Error from gpvert.m:           !!!') ;
      disp('!!!       Number of cycles < 5       !!!') ;
      disp('!!!       No movie made              !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    end ;

  end;

  else ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from gpvert.m:           !!!') ;
    disp('!!!       Number of input lags       !!!') ;
    disp('!!!       must be 1, No movie made   !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

 end ;




