c
c              INVAR2  MAIN  Program for Ward Johnson Problem
c 
c*******************************************************************************
c  This version of INVAR (completed October 1982, but modified by Bert W. Rust
c  in April 1996) is intended for use with the DISSPLA plotting package. 
c  It includes all of the original 18 subroutines (with minor adjustments) 
c  used by Rall and Funderlic (May 1980) except for PLOT.  The latter has been 
c  replaced by the subroutines PLOT1 and PLOT2 which contain calls to DISSPLA 
c  subroutines.
c-------------------------------------------------------------------------------
c
c   The following program and set of subroutines provide an interactive program
c       for solving either linear or nonlinear least squares problems.
c
c   The package consists of a MAIN program and 30 subroutines including:  
c       CH, PR, VARPRO, ORFAC1, ORFAC2, DPA, BACSUB, POSTPR, COV, VARERR, ADA,
c       STAT, ARRAYD, PLOT1 and PLOT2.  
c
c   The user is asked to supply 5 subroutines: CALC, DERIV, WT, TRANS, BAKTRN.
c
c   The method used for solving the least squares problem is described in
c       the comments preceding subroutine VARPRO.  Also included are a
c       description of the parameters, comments on convergence failures and
c       a list of references.
c
c           For further information, one should refer to the report:
c       "Interactive VARPRO (INVAR), a Nonlinear Least Squares Program", 
c           by J.E. Rall and R.E. Funderlic, ORNL/CSD-55 (May 1980).
c
c*******************************************************************************
c
      PARAMETER ( mmax = 3000 ,  miv = 6 ,  maxn = 100 ,  maxkg = 50 ,  
c
     1            maxp = 200 ,    nf = 8192 )
c                                \_________/ 
c                                     ^
c                                     |_______ must be a power of 2 .ge.  mmax
c
c                 mmax  = max. no. observations (data points)
c                 miv   = max. no. independent variables
c                 maxn  = max. no. linear parameters
c                 maxkg = max. no. nonlinear parameters
c                 maxp  = max. no. nonzero partial derivatives
c                 nf    = no. frequencies in Fourier spectrum of residuals
c                         ( must be a power of 2 .ge.  mmax )
c
c-------------------------------------------------------------------------------
c
      dimension   a(mmax,maxn+maxp+2),       dph(maxn+1,maxkg), 
     1           cr(maxn+maxkg,maxn+maxkg)
c
      dimension  t(mmax,miv), y(mmax),     w(mmax),     weight(mmax),
     1           alf(maxkg),  c(maxn),     ph(maxn+1),  eta(mmax),
     2           r(mmax),     by(mmax),    bketa(mmax), bt(mmax,miv),
     3           tempx(mmax), tempy(mmax), xpl(mmax),   ypl(mmax)
c
      dimension  freq(nf+1), perdg(nf+1), amp(nf+1),  cumper(nf+1),
     1           flb(2),     clb(2),      fub(2),     cub(2),
     2           spwork(10*nf + 9)
c
      real ndi(maxp)
c
      complex  cmplxts(mmax/2), cwrk1(4*nf), cwrk2(5*mmax)
c
      integer  p, ia(maxn+1), inc(maxkg,maxn+1)
c
      character*1 nam(40)
c
      character*4  corlab(maxn+maxkg)
c
      character*6  clab(maxn), alflab(maxkg)
c
      character*12  wavlab
c
      character*47  tlabel(miv), ylabel
c
      character*56  title
c
      common /titlab/ title
      common / tlab / tlabel
      common / ylab / ylabel
c
      character*47  btlabl(miv), bylabl
      character*56  btitle
c
      common /btilab/ btitle
      common /btlab / btlabl
      common /bylab / bylabl
c
      character*1 ians
      character*4 icode
      character*26 fname, modfnam
c
      logical  ip, iq, nc, iwt, wrnbef, ts, itrn, itp
      logical  df, ex, pmsig, whn, modfile
c
      common /blokip/ ip
      common /blkndi/ ndi
      common /blkinc/ inc
      common /ncblk/ nc
      common /blkph/ ph
      common /blkdph/ dph
      common /blkts/ ts
      common /blkitp/ itp
      common /blockp/ p
c
      common /mtotbl/ mtot
c                        --- Transfers total no. data pts. to CALC  and  DERIV
c
      common /modelbl/ ncon, ncwav, npexp, npsin
c
c                   ncon = no. constant (baseline correction) terms in model
c                  ncwav = no. coupled sine/cosine waves in model
c                  npexp = no. pure exponentials (0 or 1)
c                  npsin = no. pure (undamped) sinusoids
c
      data  itmax /50/ , eps1 /1.0e-5/ , iprint /1/
c
      data  iq, iwt, wrnbef, itrn / 4 * .false. /
c
c--------------- initialization segment of program ---------------------
c
c-----      ip = .false.
c
      ip = .true.
c
      itp = .true.
c
      nc = .false.
c
      ts = .false.
c
170   format (a1)
c
      open ( unit=10 , file = 'lpout' , status='new' )
c
      print *, ' Hard copy output will be in file  lpout  '
c
c--------------
  1   continue
c--------------
c
      print *, ' Enter file name for standard data file  '
      print *, ' (Type "null" if data not in a standard data file)  '
      read 660, fname
  660 format ( a26 )
c
      if ( fname .ne. 'null' ) then
c                                  -- Read data from standard file
             df = .true.
c
             inquire ( file = fname , exist = ex )
c
             if ( .not. ex ) then
                               print 305, fname
                     go to 1
             end if
c
             open ( unit = 4 , file = fname , status = 'old' ) 
c
             read (4,*) m
c
             read (4,20) title
             btitle = title
c
             read (4,4) tlabel(1), iv
    4        format ( a47, i3 )
             btlabl(1) = tlabel(1)
c
             if ( iv .le. 1 ) then
                               iv = 1
             else
                 do 3 il=2,iv
                    read (4,90) tlabel(il)
                    btlabl(il) = tlabel(il) 
    3            continue
             end if
c
             read (4,90) ylabel
             bylabl = ylabel
c                                                              ( y_2 )
c---   The y-array will be stored internally in the order  y = ( --- )
c                                                              ( y_1 )
c
c  where:  First array in data file = y_1 ~= Sum[cosine functions]
c
c         Second array in data file = y_2 ~= Sum[sine functions]
c
             mo2 = m / 2
             mo2p1 = mo2 + 1
c
             do 2 i = 1, mo2
                read (4,*) ( t(mo2+i,j), j=1,iv ), y(mo2+i)
    2        continue
c                     ---- reads y_1 into lower half of y-array
c
  183        print *, ' Enter gain correction factor for y_2 data '
             print *, '          (eg,  gcorr2 = 0.9894 )          '
             read *, gcorr2
c
             print 185, gcorr2
  185        format(' '/' '/' Gain correction factor for y2 =', 0pf9.6 )
             print *, ' Is this value OK? '
             read 170, ians
             if ( (ians .ne. 'y') .and. (ians .ne. 'Y') ) go to 183
             if (ip) write (10,185) gcorr2
c
             do  i = 1, mo2
                read (4,*) t(i,1), y(i)
                y(i) = gcorr2 * y(i)
             end do
c                     ---- reads y_2 into upper half of y-array
c                      --- and applies the gain correction factor
c
             rewind 4
             close ( unit = 4 , status = 'keep' )
c
      else
c           ------------- Type in title and variable names and read
c                         data from a non-standard file.
c
  190       print 200
  200       format (' no. of data points:  ') 
            read (5,*) m
            if ( m.lt.1 .or. m.gt.mmax ) go to 190 
            if (ip) write (10,210) m
  210       format (' '/' ','no. of data points: ',i3)
c
            if (iq) go to 770
c
  280       print 290
  290       format (' no. of independent variables:  ') 
            read (5,*) iv 
            if ( iv.lt.1 .or. iv.gt.miv ) go to 280
            if (ip) write (10,300) iv
  300       format (' no. of independent variables: ',i2)
c
            if (iq) go to 770
c
            print 10
   10       format (/' enter title (maximum of 56 characters)'/)
            read 20, title
   20       format (a56)
            btitle = title
c 
            print 80
   80       format (/' enter label for independent variable 1 ',
     1               '(max. of 47 characters)'/)
            read 90, tlabel(1)
   90       format (a47)
            btlabl(1) = tlabel(1) 
c
            print 100
  100       format (/' enter label for dependent variable ',
     1               '(max. of 47 characters)'/)
            read 90, ylabel
            bylabl = ylabel
c
            if (ip) write (10,1120) title
            print 180
 180        format (' are your data, weights, initial values, etc. '/
     1         '                already in a VARPRO-defined file ?  ') 
            read 170, ians
            if (ians.eq.'y') go to 640
c
c          ----------------------------------------------------------- 
c           this section allows data to be read from other files
c          -----------------------------------------------------------
c
  301       df=.false.
            print 302
  302       format(' are data in a file ?  ')
            read 170,ians 
            if(ians.ne.'y') go to 309
            df=.true.
  303       print 304
  304       format(' what is the filename ? -- no more than 26 chars'/
     *          '                           (type "quit" to exit)  ' )
            read 660, fname
            if (fname.eq.'quit') go to 301
            inquire (file=fname,exist=ex)
            if (ex) go to 306
            print 305,fname
  305       format(' ',a26,' does not exist.')
            go to 303 
  306       open(unit=4,file=fname,status='old')
c
            do 307 i=1,m
               read(4,*) (t(i,j),j=1,iv),y(i)
  307       continue
            rewind 4
            close(unit=4,status='keep')
c 
      end if
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      print *, ' Would you like to average data before analyzing it? '
      read 170, ians
c
      if ( (ians .eq. 'Y') .or. (ians .eq. 'y') ) then
c
c                                                                      ( y_2 )
c---         Recall the y-array is stored internally in the order  y = ( --- )
c                                                                      ( y_1 )
c
c            where:  First array in data file = y_1 ~= Sum[cosine functions]
c
c                    Second array in data file = y_2 ~= Sum[sine functions]
c
             print *, ' Enter number of points to average. '
             read (5,*) nav
c
             morig = m
             mo2or = m / 2
c
             m = m / nav
             if ( (m / 2)*2 .ne. m ) m = m - 1
             mo2 = m / 2
c
             if (ip) write (10,222) nav, m
c
             do i = 1, mo2
c
                t2av = 0.0
                y2av = 0.0
c
                do k = 1, nav
c 
                   ifr2 = (i - 1)*nav + k
                   t2av = t2av + t(ifr2,1)
                   y2av = y2av + y(ifr2)
c
                end do
c
                t(i,1) = t2av / float( nav )
                y(i)   = y2av / float( nav )
c
                if (ip) write (10,223) i, t(i,1), y(i)
c
             end do
c
c
             if (ip) write (10,224) nav, m 
c
             do i = 1, mo2
c
                t1av = 0.0
                y1av = 0.0
c
                do k = 1, nav
c
                   ifr1 = (mo2or) + (i - 1)*nav + k
                   t1av = t1av + t(ifr1,1)
                   y1av = y1av + y(ifr1)
c
                end do
c
                t(mo2+i,1) = t1av / float( nav )
                y(mo2+i)   = y1av / float( nav )
c
                if (ip) write (10,223) i, t(mo2+i,1), y(mo2+i)
c
             end do
c
      end if
c
  222 format('1'/' '/' Using ', i2, '-point averaging, new  m =', i4,
     1         / ' '/'   i', 8x, 't(i)', 11x, 'y_2(i)' / ' ' )
c
  223 format(' ', i4, 1p2e15.6 )
c
  224 format('1'/' '/' Using ', i2, '-point averaging, new  m =', i4,
     1         / ' '/'   i', 8x, 't(i)', 11x, 'y_1(i)' / ' ' )
c
c--- Plot the Input Data
c
      mo2 = m / 2
      mo2p1 = mo2 + 1
      do i = 1, m
         tempy(i) = 0.0
      end do
      ts = .true.
c   -------------------------------------------------------------------
      call TWOPLT ( t(1,1), y, mo2, 
     1              t(1,1), tempy, mo2, 0, 1,
     2              'Sines, '//title, tlabel(1), ylabel,
     3              t(mo2p1,1), y(mo2p1), mo2,
     4              t(mo2p1,1), tempy(mo2p1), mo2, 0, 1,
     5              'Cosines, '//title, tlabel(1), ylabel, ts )
c   -------------------------------------------------------------------
c
c
      print *, ' Compute real periodograms of the two wave-forms? '
      read 170, ians
c                                                  +++++++++++
      if ( (ians .ne. 'y') .and. (ians .ne. 'Y') )  go to 227
c                                                  +++++++++++
c
c-- Calculate periodograms of the two real waves separately
c
      do iperd = 1, 2
c
c
         if ( iperd .eq. 1 ) then
                                  ip1 = mo2 + 1
                                  wavlab = 'Cosine Wave '
                  else
                       ip1 = 1
                       wavlab = ' Sine Wave  '
         end if
c
         mo2 = m / 2
c
         delt =  ( t(mo2,1) - t(1,1) ) / float(mo2 - 1)
c
c
         if ( iperd .eq. 1 ) then
c
               do i = 1, mo2
                  cmplxts(i) = CMPLX( y(mo2+i) , 0.0 )
               end do               
c
c            ----------------------------------------------------------
               call FFTPLOT ( t, cmplxts, mo2, delt, nf, freq, cwrk1,
c
     1                        cwrk1(nf+1), cwrk1(2*nf+1), spwork, 
c
     2                        spwork(nf+1), 'Cosine Wave Form', 
c
     3                        spwork(2*nf+1), spwork(3*nf+1),
c
     4                        spwork(4*nf+1), spwork(5*nf+1), ts )
c            -----------------------------------------------------------
c            on return:  spwork(j), j = 1   , nf,   contains real component,
c                      & spwork(j), j = nf+1, 2*nf, contains imag. component
c                                                   of DFT[ Cosine wave form ]
         end if
c
c
         if ( iperd .eq. 2 ) then
c
               do i = 1, mo2
                  cmplxts(i) = CMPLX( 0.0 , y(i) )
               end do               
c
c            ----------------------------------------------------------
               call FFTPLOT ( t, cmplxts, mo2, delt, nf, freq, cwrk1,
c
     1                        cwrk1(nf+1), cwrk1(2*nf+1), spwork, 
c
     2                        spwork(nf+1), 'i*(Sine Wave Form)', 
c
     3                        spwork(2*nf+1), spwork(3*nf+1),
c
     4                        spwork(4*nf+1), spwork(5*nf+1), ts )
c            -----------------------------------------------------------
c            on return:  spwork(j), j = 1   , nf,   contains real component,
c                      & spwork(j), j = nf+1, 2*nf, contains imag. component
c                                                   of DFT[ i*(Sine wave form) ]
         end if
c
c
c      ------------------------------------------------------------------
         call PERDGR ( y(ip1), mo2, delt, nf, spwork, freq, perdg, amp,
c
     1                 sump, pave, pmax, fpmax, fisher, pmsig,
c
     2                 cumper, nout, flb, clb, fub, cub, whn, info )
c      ------------------------------------------------------------------
c
         if ( info .gt. 0 ) then
c
                print *, ' Error in call to PERDGR,  info =', info,
     1                   '     iperd =', iperd
                if (ip) write (10, 1980) info, iperd
c
c              +++++++++++
                go to 227
c              +++++++++++                
c
         end if
c
         xpl(1) = flb(1)
         ypl(1) = clb(1)
         xpl(2) = flb(2)
         ypl(2) = clb(2)
         xpl(3) = flb(2)
         ypl(3) = 1.0
         xpl(4) = fub(2)
         ypl(4) = cub(2)
         xpl(5) = fub(1)
         ypl(5) = cub(1)
c
         permx = 1.0 / fpmax
         ratio = pmax / pave
c
         if (ip) write (10,1935) wavlab, sump, pave
                     print 1935, wavlab, sump, pave
         if (ip) write (10,1945) fpmax, permx, pmax
                     print 1945, fpmax, permx, pmax
         if (ip) write (10,1940) ratio, fisher
                     print 1940, ratio, fisher
         if ( pmsig ) then 
                           write (10, 1942)
                           print 1942
              else
                   write (10,1943)
                   print 1943
         end if
c
         if ( whn .eq. .true. ) then
                                    if (ip) write (10, 1985) nout, nf
                                                 print 1985, nout, nf
                 else 
                     write (10, 1986) nout, nf
                          print 1986, nout, nf        
         end if
c
         do i = 1, nf
            spwork(i) = 0.0
         end do
c
c      -----------------------------------------------------------------
         call TWOPLT ( freq, perdg, nf, freq, spwork, nf, 0, 1, 
     1                 wavlab//'Data Periodogram', 
     1                 'Frequency', 'Variance', 
     2                 freq, cumper, nf, xpl, ypl, 5, 0, 1,
     3                 wavlab//'Cumulative Periodogram', 
     4                 'Frequency', 'Cumulative Fraction', ts )
c      -----------------------------------------------------------------
         call TWOPLT ( freq, perdg, nf, freq, spwork, nf, 0, 1, 
     1                 wavlab//'Data Periodogram', 
     1                 'Frequency', 'Variance', 
     2                 freq, amp, nf, freq, spwork, nf, 0, 1,
     3                 wavlab//'Data Amplitude Spectrum', 
     4                 'Frequency', ylabel, ts )
c      ----------------------------------------------------------------
c
c     
      end do
c
c
c++++++++++++++
  227 continue
c++++++++++++++
c
c
c
c-- Calculate periodogram of complex data obtained by combining the two waves --
c
      mo2 = m / 2
c
      delt =  ( t(mo2,1) - t(1,1) ) / float(mo2 - 1)
c
      do i = 1, mo2
c----         cmplxts(i) = CMPLX( y(i), y(mo2+i) )
         cmplxts(i) = CMPLX( y(mo2+i), y(i) )
      end do
c
      iou = 10
c
      print *, ' Plot real and imaginary components '
      print *, '                         of complex time series? '
      read 170, ians
c
      if ( (ians .eq. 'Y') .or. (ians .eq. 'y') ) then
c
c            ----------------------------------------------------------
              call FFTPLOT ( t, cmplxts, mo2, delt, nf, freq, cwrk1,
c
     1                       cwrk1(nf+1), cwrk1(2*nf+1), spwork, 
c
     2                       spwork(nf+1), 'Complex Raw Data', 
c
     3                       spwork(2*nf+1), spwork(3*nf+1),
c
     4                       spwork(4*nf+1), spwork(5*nf+1), ts )
c            -----------------------------------------------------------
c            on return:  spwork(j), j = 1   , nf,   contains real component,
c                      & spwork(j), j = nf+1, 2*nf, contains imag. component
c                                                   of the complex DFT
c
      end if
c            
c    --------------------------------------------------------------------
      call CMPLXPER ( t, cmplxts, mo2, delt, nf, freq, perdg, amp,
c
     1                spwork, spwork(nf+1), ylabel, ts,
c
     1                cwrk1, cwrk1(nf+1), cwrk1(2*nf+1), cwrk1(3*nf+1),         
c
     2                cwrk2, cwrk2(mmax+1), cwrk2(2*mmax+1),
c
     3                      cwrk2(3*mmax+1), cwrk2(4*mmax+1), iou,
c
     4                     'Raw Data ' )
c    --------------------------------------------------------------------
c
c
      ip=.true.
c
      if (ip) write (10,1120) title
c
c++++++++++++++
  220 continue
c++++++++++++++
c
      print *, ' Read model specs. and initial estimates from a file? '
      read 170, ians
      if ( (ians .eq. 'y') .or. (ians .eq. 'Y') ) then
                                                      modfile = .true.
             print *, ' Enter file name: '
             read *, modfnam
             open ( unit = 16, file = modfnam, status = 'old' )
      end if
c
      if ( .not. modfile ) then
             print *, ' Do you want baseline-correction const. terms? '
             print *, '    -- i.e., (B = 0 and C = 0) terms? '
             read 170, ians
        else
             read (16,*) ians
      end if
c           
      if ( (ians .eq. 'y') .or. (ians .eq. 'Y') ) then
                                                  ncon = 2
                                                  n    = 2
                               if (ip) write (10,228)
                               if (ip) write (10,240) n
              else
                   ncon = 0
                   n    = 0
      end if
c
  228 format(' '/' Model contains constant terms in both waveforms.')
c
      if ( .not. modfile ) then      
                               print 230
                               read (5,*) ncwav
         else
              read (16,*) ncwav
      end if
c
  230 format (' Enter no. coupled sine/cosine waves (.le. 5):  ')      
c
      n = ncon + 2 * ncwav
c                                          +++++++++++
      if ( (n .lt. 0) .or. (n .gt. maxn) )  go to 220 
c                                          +++++++++++
      if (ip) write (10,238) ncwav
  238 format (' '/' Model contains ', i2, ' coupled sine/cosine waves' )
c
      if (ip) write (10,240) n
240   format (' '/' No. of linear parameters: ',i2)
c
      kg = 2 * ncwav
c                                             +++++++++++
      if ( (kg .lt. 0) .or. (kg .gt. maxkg) )  go to 220
c                                             +++++++++++
      if (ip) write (10,270) kg
270   format (' No. of nonlinear parameters: ',i2)
c
      if ( .not. modfile ) then
               print *, ' Do you want pure exponential (C = 0) terms?  '
               print *, '  i.e. Is the reference frequency in signal?  '
               read 170, ians
         else
              read (16,*) ians
      end if
c
      if ( (ians .eq. 'y') .or. (ians .eq. 'Y') ) then
                                               npexp = 1
                                               n =  n  + 2
                                               kg = kg + 1
                               if (ip) write (10,275)
                               if (ip) write (10,240) n
                               if (ip) write (10,270) kg
          else 
               npexp = 0
               if (ip) write (10,*)  ' '
               if (ip) write (10,*) ' Reference signal not in the data',
     1                              ' (No pure exponential in model).'
c
      end if
c
  275 format(' '/' Model contains a pure exponential, i.e., the ',
     1        'reference signal is in the data.')
c
      if ( .not. modfile ) then
               print *, ' Enter no. of pure (undamped) sinusoids:  '
               read (5,*) npsin
         else
               read (16,*) npsin
      end if
c
      n  = n + 2 * npsin
      kg = kg + npsin
c
      if ( npsin .gt. 0 ) then
               if (ip) write (10,277) npsin
               if (ip) write (10,240) n
               if (ip) write (10,270) kg
          else
              npsin = 0
              if (ip) write (10,*)  ' '
              if (ip) write (10,*) ' No pure sine/cosine waves in model'
      end if
c
  277 format (' '/' Model contains ', i2, ' pure (undamped) ',
     1            'sine/cosine waves' )
c
      if (iq) go to 770
c
      go to 370
c
c-------------------------------------------------------------------
c
309   continue

      print 310
310   format (' Are weights to be entered with data ?  ')
      read 170, ians
      if (ians.ne.'y') go to 370
      print 320
      if (ip) write (10,320)
320   format (/' data points:  list independent variables, dependent',
     2        /' variable & weight = w squared, separated by blanks', 
     3        ' or commas'/)
      do 360 j=1,m
         print 330, j
330      format (' ',i3,'  ') 
         read (5,*) (t(j,i),i=1,iv),y(j),weight(j)
340      format (e14.8)
         if (ip) write (10,350) j,(t(j,i),i=1,iv),y(j),weight(j)
350      format (' ',i3,3x,1p4e14.6,/6x,1p4e14.6) 
360   continue
      go to 460
c
c-------------------------------------------------------------------------------
c
 370  continue
c
c----      print 380
c----380   format (' are weights to be calculated in subroutine WT ?  ')
c----      read 170, ians
c----      if ( ians .eq. 'y' ) go to 410
c
      do 390 j=1,m
390      weight(j)=1.0
      print 400
      if (ip) write (10,400)
400   format (' '/' Weights are set to 1  '/)
c
      go to 430
c
410   continue
c
      iwt=.true.
      if (ip) write (10,420)
420   format (' Weights calculated in subroutine WT '/)
c
430   continue
      if (df) go to 798
c
      print 440
440   format (/' data points:  list independent variables and'/
     2        ' dependent variable, separated by blanks or commas'/)
      do 450 i=1,m
         print 330, i
         read (5,*) (t(i,j),j=1,iv),y(i)
         if (ip) write (10,350) i,(t(i,j),j=1,iv),y(i)
450   continue
c
798   continue
c
c----      print 800
c----800   format (' '/' Do dependent or independent variables'/
c----     2        ' need to be transformed before analysis ?  ')
c----      read 170, ians
c----      if ( ians .ne. 'y' ) then
c
           btitle = title
           do il = 1, iv
              btlabl(il) = tlabel(il)
           end do
           bylabl = ylabel 
c    
           go to 460
c
c----      end if
c
451   itrn=.true.
c
c    ------------------------------------------
      call TRANS ( t, bt, y, by, iv, m, mmax )
c    ------------------------------------------
c
      if (.not.iq) goto 460
      print 459
459   format (/' if weights were calculated in subroutine WT,'/
     1        ' they must be recalculated with transformed data.')
c
      go to 780
c
c              --------------------------
460   if (iwt)  call WT ( m, weight, y )
c              --------------------------
      if ( kg .eq. 0 ) then
                       iq = .true.
             go to 770 
      end if
c
      print 470
      if (ip) write (10,470)
470   format (' '/' Enter initial values of nonlinear parameters:')
c
      do kwav = 1, ncwav
         kalf = 2*kwav - 1
         if ( .not. modfile ) then
                              print 480, (kwav + ncon/2)
                              read (5,*) alf(kalf)
                              print 481, (kwav + ncon/2)
                              read (5,*) alf(kalf+1)
                 else
                     read (16,*) alf(kalf), alf(kalf+1)
         end if
      end do
  480 format(' Enter estimate of B', i1, ':  ')
  481 format(' Enter estimate of C', i1, ':  ')
c
      if (npexp .gt. 0) then
              kwav = (ncon/2) + ncwav + 1
              kalf = 2*ncwav + 1
              if ( .not. modfile ) then
                                       print 482, kwav
                                       read (5,*) alf(kalf)
                     else
                          read (16,*) alf(kalf)
              end if
      end if
  482 format(' Enter estimate of B', i1, 
     1               ' for pure exponential:  ' )
c
      if (npsin .gt. 0) then
                            do isin = 1, npsin
                               kwav = (ncon/2) + ncwav + npexp + isin
                               kalf = 2*ncwav + npexp + isin
                               print 483, kwav
                               read (5,*) alf(kalf)
                            end do
      end if
c
  483 format(' Enter estimate of C', i1, ' for undamped sine/cosine: ')
c
c    --------------------------------------------------------------
      call PR ( 99999, 'palf', mmax, m, n, kg, p, iv, t, y, weight,
c
     *                                                       alf ) 
c    --------------------------------------------------------------
c
c--------------
 491  continue
c--------------
c
c  determine incidence matrix 
c
      np1=n+1
      ph(np1)=0.
      do 500 j=1,np1
         do 500 k=1,kg
            inc(k,j)=0
  500 continue
c
      do kwav = 1, ncwav
         irow = 2*kwav - 1
         jcol = ncon + 2*kwav - 1
         inc(irow,jcol)     = 1 
         inc(irow,jcol+1)   = 1
         inc(irow+1,jcol)   = 1
         inc(irow+1,jcol+1) = 1
      end do
c
      if (npexp .gt. 0) then
                            irow = 2*ncwav + 1
                            jcol = ncon + 2*ncwav + 1
            inc(irow,jcol)   = 1
            inc(irow,jcol+1) = 1
      end if
c
      if (npsin .gt. 0) then
                            do isin = 1, npsin
                               irow = 2*ncwav +   npexp + isin
                               jcol = ncon + 2*ncwav + 2*npexp 
     1                                               + (2*isin - 1)
                            end do
            inc(irow,jcol)   = 1
            inc(irow,jcol+1) = 1
      end if
c
c    --------------------------------------------------------------
      call PR ( 99999, 'pinc', mmax, m, n, kg, p, iv, t, y, weight,
c
     *                                                       alf ) 
c    --------------------------------------------------------------
c
      p=0 
c
      do 540 j=1,np1
         do 540 i=1,kg
540         if (inc(i,j).eq.1) p = p + 1
c
      if (iq) go to 770
c
      iq = .true.
c
c  initialize numerical derivative increments
c
c----      print 550
c----550   format (' all derivatives in closed form ?  ')
c----      read 170, ians
      ians = 'y'
c
      if (ians.eq.'y') go to 780
c
      print 560
560   format (' all numerical derivative increments the',' same ?  ')
      read 170, ians
      if (ians.eq.'y') go to 600
      print 570
570   format (' numerical derivative increments'/ 
     2        ' (enter 0 for those to be found from closed form)')
      if (ip) write (10,580)
580   format (' '/' ','numerical derivative increments')
      do 590 k=1,p
         print 330, k
         read (5,*) ndi(k)
590      if (ip) write (10,480) k,ndi(k)
      go to 780
600   continue
      print 610
610   format (' increment values:  ')
      read (5,*) ans
      if (ip) write (10,620) ans
620   format (' ',' ','increment values:  ',e15.7)
      do 630 i=1,p
630      ndi(i)=ans 
      go to 780
c
c--------------
 640  continue
c--------------
c
c  read in data that has been saved
c
      print 650
650   format (' type name of file (max. of 26 characters):  ')
659   read 660, fname
      inquire (file=fname,exist=ex)
      if (ex) goto 662
      print 661,fname
661   format (/' ',a26,' does not exist.  try again.')
      goto 659
662   open(unit=3,file=fname,status='old')
      read (3,670) (nam(i),i=1,11),m,(nam(i),i=1,20),n,(nam(i),i=1,23),
     2             kg,(nam(i),i=1,24),iv,(nam(i),i=1,32)
670   format (11a1,i4/20a1,i4/23a1,i4/24a1,i4/32a1)
      do 680 i=1,m
         read (3,340) (t(i,j),j=1,iv),y(i),weight(i)
 680  continue
c
      read (3,690) (nam(i),i=1,23),(alf(i),i=1,kg)
690   format ( 23a1 / 10e14.8 )
c
      np1=n+1
      ph(np1)=0.
      do 700 j=1,np1
         do 700 k=1,kg
700         inc(k,j)=0
      do 720 j=1,kg 
         read (3,710) (nam(i),i=1,39),(ia(k),k=1,np1)
710      format (39a1,11i2)
         do 720 i=1,np1
            if (ia(i).eq.0) go to 720
            ib=ia(i)
            inc(j,ib)=1
720   continue
      p=0 
      do 730 j=1,np1
         do 730 i=1,kg
730         if (inc(i,j).eq.1) p=p+1
      iq=.true.
      read (3,740) (nam(i),i=1,34),(ndi(i),i=1,p) 
740   format (34a1 / 10e14.8)
      rewind 3
      close(unit=3,status='keep')
      print 750
750   format (/' would you like your initial',
     2        ' values file printed ?  ')
      read 170, ians
      if (ians.ne.'y') go to 780
c
c    ----------------------------------------------------------------
      call PR ( 1, icode, mmax, m, n, kg, p, iv, t, y, weight, alf )
c    ----------------------------------------------------------------
c
c--------------
 780  continue
c--------------
c 
      print *, ' '
c
c
c----------------------------- Command point -----------------------------------
c
c    ++++++++++
770   continue
c    ++++++++++
c
      if (wrnbef) go to 840
c
c     the following questions are asked only once per execution
c
c++++++++++++++
 790  continue
c++++++++++++++
c
c-----790   print 820
c-----820   format (/' does data follow an evenly spaced time series ?  ')
c-----      read 170, ians
c-----      if ( ians .eq. 'y' ) ts = .true.
c
      ts = .true.
c
      print 830
830   format (/' '/'   note:   the command "sc" will save',
     2        ' your data file at any time')
c
c
      wrnbef=.true. 
c
840   continue
c
      print 850
850   format (' '/' next step ? (type "h" for help)  ')
      read 860, icode
860   format (a4)
c
      if (icode.ne.'h') go to 950
c
872   print 870
      if (ip) write (10,870)
870   format ('1'/' the command codes fall into four main categories'/
     2  /'     print-  to examine current values of initialized data'/
     3   '     change- to correct any initialized value, transform'/
     4   '             data, and either add or delete data points'/
     5   '     run-    to run the least squares program to find the'/
     6   '             optimal parameters, or to run a statistical'/
     7   '             analysis of the data returned'/
     8   '     save/   to save current data or '/
     9   '     stop-   to stop the program, with or without'/ 
     *   '             saving the initialized data.'//
     1   ' type "p","c","r", or "s" for the specific commands in'/
     2   ' the above categories.   type "a"',' to obtain a listing of'/         
     3   ' all command codes.   ')
c
      read 170, ians
c
      if ( ians.eq.'c' ) go to 890
c
      if ( ians.eq.'r' ) go to 910
c
      if ( ians.eq.'s' ) go to 930
c
      if ( ians.ne.'p' .and. ians.ne.'c' .and. ians.ne.'r' 
     *                 .and. ians.ne.'s' .and. ians.ne.'a' ) go to 872       
c
      print 880
      if (ip) write (10,880)
880   format (/' '/' command codes for printing'//
     2   '     pm    print the no. of data points'/ 
     3   '     pn    print the no. of linear parameters'/
     4   '     pkg   print the no. of nonlinear parameters'/
     5   '     piv   print the no. of independent variables'/ 
     6   '     ptw   print the data points & weights'/
     7   '     palf  print the values of the nonlinear',' parameters'/
     8   '     pndi  print the numerical derivative',' increments'/
     9   '     pinc  print which terms depend on nonlin',' parameters'/ 
     *   '     lp    turn on the line printer flag'/
     1   '     np    turn off the line printer flag'/
     2   '     tp    turn on the flag to plot on crt'/
     3   '     ntp   turn off the flag to plot on crt'//)
c
      if (ians.ne.'a') go to 770
c
890   continue
      print 900
      if (ip) write (10,900)
900   format (/' '/' command codes for changing'//
     2    '     cm    change no. of data points'/
     3    '     cn    change no. of linear parameters'/
     4    '     ckg   change no. of nonlinear parameters'/
     5    '     civ   change no. of independent variables'/
     6    '     ct    change one or more data points'/
     7    '     cw    change one or more weights'/
     8    '     calf  change one or more nonlinear',' parameters'/
     9    '     cndi  change one or more numer deriv',' increments'/
     *    '     ctrn  transform the data'/ 
     *    '     cinc  change which terms depend on nonlin parameters'//)
c
      if (ians.ne.'a') go to 770
c
910   continue
      print 920
      if (ip) write (10,920)
920   format (/' '/' command codes for running'// 
     2    '     rv    run the least squares program'/
     3    '     rs    run a statistical analysis of',' results'/
     4    '     fsr   run Fourier spectral analysis of residuals'//)
c
      if (ians.ne.'a') go to 770
c
930   continue
      print 940
      if (ip) write (10,940)
940   format (/' '/' command codes to save data file and/or stop'//
     2    '     sc    save current values of data & continue '/
     3    '     sd    save current values of data & stop'/
     4    '     sr    save current values of indep. variables ',
     5    'and residuals & continue'/'     s0    stop without saving'//)
c
      go to 770
c
c-------------------------------------------------------------------------------
c
950   continue
c
      if ( icode .eq. 'rv' ) go to 960
c
      if ( icode .eq. 'rs' ) go to 1580
c
      if ( icode .eq. 'fsr' ) go to 1900
c
      if ((icode.eq.'pm').or.(icode.eq.'pn').or.(icode.eq.'pkg').or.
     2    (icode.eq.'piv').or.(icode.eq.'ptw').or.(icode.eq.'palf').or.
     3    (icode.eq.'pndi').or.(icode.eq.'lp').or.(icode.eq.'np').or. 
     4    (icode.eq.'pinc').or.(icode.eq.'tp').or.(icode.eq.'ntp'))
     5                                                       go to 1600 
c
      if ( icode .eq. 'cm' )  then
  945                             print 200
                                  read (5,*) m
                     if ( m.lt.1 .or. m.gt.mmax ) go to 945 
      end if
c
      if ( icode .eq. 'cn' )  then
  947                            print 230
                                 read (5,*) n
                                 if ( n.lt.0 .or. n.gt.maxn ) go to 947 
                                 if (ip) write (10,240) n
      end if
c
      if ( icode .eq. 'ckg' )  then
  948                             print 260
                                  read (5,*) kg 
                              if ( kg.lt.0 .or. kg.gt.maxkg ) go to 948
                                  if (ip) write (10,270) kg
      end if
c
260   format (' no. of nonlinear parameters:  ')
c
      if ( icode .eq. 'civ' ) then
  949                              print 290
                                   read (5,*) iv
                               if ( iv.lt.1 .or. iv.gt.miv ) go to 949 
                                   if (ip) write (10,300) iv
      end if
c
      if ( icode .eq. 'ctrn' ) go to 451
c
      if ( icode .eq. 'cinc' ) go to 491
c
      if ( (icode.eq.'ct') .or. (icode.eq.'cw') .or. (icode.eq.'calf') 
     2                     .or. (icode.eq.'cndi') )
     *                                              go to 1590 
c
      if ( icode .eq. 'sc')  go to 1680
c
      if ( icode .eq. 'sd' ) go to 1680
c
      if ( icode .eq. 'sr' ) go to 1630
c
      if ( icode .eq. 's0' ) go to 1790
c
      go to 1610
c
c
c------------------------- Run the least squares program -----------------------
c
960   continue
c
      if ( kg .eq. 0) go to 1090 
c
      iprint = 1
970   print 980
980   format (/' do you wish to change any defaults ',
     2        'in the least squares program ?'/
     3        '        (     50 : max. no. iterations   )'/
     4        '        ( 1.0e-5 : converg. tolerance    )'/
     5        '        (      1 : output printing index )   ')
      read 170, ians
      if (ians.ne.'y') go to 1110
c
990   print 1000
1000  format (/' command codes for changing defaults'/
     2        '     cmni  change max. no. iterations'/
     3        '     cct   change converg. tolerance'/
     4        '     copi  change output printing index  ')
      read 860, icode
c
      if ( icode .eq. 'cmni') go to 1010
c
      if ( icode .eq. 'cct' ) go to 1050
c
      if ( icode .eq. 'copi' ) go to 1070
c
      go to 970 
c
1010  print 1020
1020  format (/' the number of iterations must be greater '/
     2         '     than 0 and less than or equal to 250.'/
     3         '                       how many iterations ?   ')
      read (5,*) itmax
      if ( itmax .le. 0  .or. itmax .gt. 250)  go to 1010
1030  print 1040
1040  format (/' do you wish to change any others ?  ')
      read 170, ians
      if ( ians .ne. 'y' ) go to 1110
      go to 990
c
1050  print 1060
1060  format (/' the convergence criterion must not be less ',
     2        'than 1.5e-8. '/1h ,
     3        'type in the new convergence criterion.    ') 
      read (5,*) eps1
      if ( eps1 .lt. 1.5e-8 ) go to 1050
      go to 1030
c
1070  print 1080
1080  format (/' type output printing index: ')
      read (5,*) iprint
1090  if ( kg .eq. 0 ) iprint = 0
      if (ip) write (10,1100) iprint
1100  format (' '/' ','output printing index:',i3)
      if ( kg .ne. 0 ) go to 1030 
c
1110  continue
      print 1120, title
      if (ip) write (10,1120) title
1120  format ('1'/' '/' '/' ', a56 /) 
c
      mtot = m
c    ------------------------------------------------------
      call VARPRO ( n, kg, m, mmax, n+p+2, iv, t, y,
c
     *              w, weight, a, iprint, alf, c,
c
     *              ierr, itmax, eps1 )
c    ------------------------------------------------------
c
      if ( ierr .le. -4 ) go to 770
      nc = .false.
      npkg = n + kg
c
      do 1180 i=1,m 
         ii=i
         mtot = m
c       ------------------------------------------------------
         call CALC ( ii, kg, n, mmax, iv, c, alf, t, eta(i) )
c       ------------------------------------------------------
 1180 continue
c
      do 1290 i=1,m 
         r(i) = y(i) - eta(i)
 1290 continue
c
c-------------------------------------------------------------------------------
c
c  Plot residuals and fit with data
c
      do 1480 inv=1,iv
c
         do 1380 i=1,m 
            tempx(i)=t(i,inv)
1380     continue
c
c
c---------  Plot residuals and fit on same VG graph 
c
         do 1355 i=1,m
            tempy(i) = 0.0
 1355    continue
c
         mo2 = m / 2
         mo2p1 = mo2 + 1
c
         if (ts.and.itp) then
c               --------------------------------------------------------
                 call SCPLT2 ( t(1,inv), r, mo2, 
     *                         t(1,inv), tempy, mo2, 0, 1,
     *                           'Sine Wave Residuals (data - model)',
     *                              tlabel(inv), ylabel )
c               --------------------------------------------------------
                 call SCPLT2 ( t(mo2p1,inv), r(mo2p1), mo2, 
     *                         t(mo2p1,inv), tempy(mo2p1), mo2, 0, 1, 
     *                         'Cosine Wave Residuals (data - model)',
     *                              tlabel(inv), ylabel )
c               --------------------------------------------------------
         end if
c
         if (.not.ts) then
c               --------------------------------------------------------
                 call SCPLT2 ( t(1,inv), r, mo2, 
     *                         t(1,inv), tempy, mo2, 17, 1, 
     *                         'Sine Wave Residuals (data - model)',
     *                              tlabel(inv), ylabel )
c               ---------------------------------------------------------
                 call SCPLT2 ( t(mo2p1,inv), r(mo2p1), mo2, 
     *                         t(mo2p1,inv), tempy(mo2p1), mo2, 17, 1, 
     *                         'Cosine Wave Residuals (data - model)',
     *                              tlabel(inv), ylabel )
c               ---------------------------------------------------------
         end if
c
         if ( iv .eq. 1 ) then
c                             -- Set up x and y arrays to plot y-fit
               tmin = tempx(1)
               tmax = tempx(1)
               do 1421 i=2,m 
               if (tempx(i) .lt. tmin) tmin = tempx(i)
               if (tempx(i) .gt. tmax) tmax = tempx(i)
 1421          continue
c
               mtot = mmax
               mmxo2 = mmax / 2
               mmxo2p1 = mmxo2 + 1
c
               delta = (tmax - tmin) / float(mmxo2 - 1) 
c
               do 1430 i = 1, mmxo2
c
                  xpl(i) = tmin + (i-1)*delta
                  xpl(mmxo2+i) = xpl(i)
c
 1430          continue
c
               do 1435 i = 1, mmax
c
c                ------------------------------------------------------
                  call CALC ( i, kg, n, mmax, 1, c, alf, xpl, ypl(i) )
c                ------------------------------------------------------
c
 1435          continue
c
               mo2 = m / 2
               mo2p1 = mo2 + 1
c             --------------------------------------------------------          
               call SCPLT2 ( t(1,inv), y, mo2, 
     1                       xpl, ypl, mmxo2, 17, 1, 
     2                       'Sines, '//title, tlabel(inv), ylabel )
c             --------------------------------------------------------
               call SCPLT2 ( t(mo2p1,inv), y(mo2p1), mo2, 
     1                       xpl(mmxo2p1), ypl(mmxo2p1), mmxo2, 17, 1, 
     2                       'Cosines, '//title, tlabel(inv), ylabel )
c             --------------------------------------------------------
c
               if ( ts .and. itp ) then
c                ------------------------------------------------------
                  call TWOPLT ( t(1,inv), y, mo2, 
     1                          xpl, ypl, mmxo2, 17, 1,  
     2                          'Sines, '//title, tlabel(inv), ylabel,
     3                          t(1,inv), r, mo2, 
     3                          t(1,inv), tempy, mo2, 0, 1, 
     4                          'Sine Wave Resids (data - model)',
     5                              tlabel(inv), ylabel, ts )
c                ------------------------------------------------------
                  call TWOPLT ( t(mo2p1,inv), y(mo2p1), mo2, 
     1                         xpl(mmxo2p1), ypl(mmxo2p1), mmxo2, 17, 1,
     2                          'Cosines, '//title, tlabel(inv), ylabel,
     3                          t(mo2p1,inv), r(mo2p1), mo2, 
     3                          t(mo2p1,inv), tempy(mo2p1), mo2, 0, 1, 
     4                          'Cosine Wave Resids (data - model)',
     5                              tlabel(inv), ylabel, ts )
c                ------------------------------------------------------
               end if
c
               if ( .not. ts ) then
c                ------------------------------------------------------
                  call TWOPLT ( t(1,inv), y, mo2, 
     1                          xpl, ypl, mmxo2, 17, 1,  
     2                          'Sines, '//title, tlabel(inv), ylabel,
     3                          t(1,inv), r, mo2, 
     3                          t(1,inv), tempy, mo2, 17, 1, 
     4                          'Sine Wave Resids (data - model)',
     5                              tlabel(inv), ylabel, ts )
c                ------------------------------------------------------
                  call TWOPLT ( t(mo2p1,inv), y(mo2p1), mo2, 
     1                         xpl(mmxo2p1), ypl(mmxo2p1), mmxo2, 17, 1,
     2                          'Cosines, '//title, tlabel(inv), ylabel,
     3                          t(mo2p1,inv), r(mo2p1), mo2, 
     3                          t(mo2p1,inv), tempy(mo2p1), mo2, 17, 1, 
     4                          'Cosine Wave Resids (data - model)',
     5                              tlabel(inv), ylabel, ts )
c                ------------------------------------------------------
               end if
c
         end if
c
1440     if ( .not. itrn ) go to 1480
c
c---     Plot back transformed data
c
c       -------------------------------
         call BAKTRN ( eta, bketa, m )
c       -------------------------------
c
         if ( iv .eq. 1 ) then
c                             ------------------------------
                               call BAKTRN ( ypl, ypl, mmax)
c                             ------------------------------
         end if
c
         do 1460 i=1,m 
            tempx(i) = bt(i,inv)
 1460    continue
c
         if ( iv .eq. 1 ) then
c                             -- Plot Back-transformed y-fit
               tmin = tempx(1)
               tmax = tempx(1)
               do 1461 i=2,m 
                  if (tempx(i) .lt. tmin) tmin = tempx(i)
                  if (tempx(i) .gt. tmax) tmax = tempx(i)
 1461          continue
c
               mmxo2 = mmax / 2
               mmxo2p1 = mmxo2 + 1
c
               delta = (tmax - tmin) / float(mmxo2 - 1)
c
               do 1470 i = 1, mmxo2
c
                  xpl(i) = tmin + (i-1)*delta
                  xpl(mmxo2+i) = xpl(i)
c
 1470          continue
c
               mo2 = m / 2
               mo2p1 = mo2 + 1
c             ---------------------------------------------------------
               call SCPLT2 ( bt, by, mo2, 
c
     1                       xpl, ypl, mmxo2, 17, 1, 
c
     1                       'Sines, '//btitle, btlabl(inv), bylabl )
c             ----------------------------------------------------------
               call SCPLT2 ( bt(mo2p1,1), by(mo2p1), mo2, 
c
     1                       xpl(mmxo2p1), ypl(mmxo2p1), mmxo2, 17, 1, 
c
     1                       'Cosines, '//btitle, btlabl(inv), bylabl )
c             ----------------------------------------------------------
c
c
         end if
c
 1480 continue
c
c-------------------------------------------------------------------------------
c
      print 1130
1130  format (' Print table of observed & predicted values?  ')
      read 170, ians
c
      if ( ians .eq. 'y' ) then
c                              -- Print table of observed and predicted values
            do 1212 inv=1,iv
               print 1190, inv
               if (ip) write (10,1190) inv 
               do 1210 i=1,m 
                  print 1200, i,t(i,inv),y(i),eta(i),a(i,npkg+1),
     1                        weight(i)
                  if (ip) write (10,1200) i,t(i,inv),y(i),eta(i),
     2                                    a(i,npkg+1), weight(i)
1210           continue
1212        continue
c
      end if
c
1190  format ('1'/' '/' '/
     1        ' index',2x,'ind var', i2, 4x,'obs val',6x,'pred val',5x,
     2        'wtd resid',5x,'weights'/)
1200  format (' ',i3,1p5e13.5)
c
c-------------------------------------------------------------------------------
c 
      if ( itrn ) then
c
            print 1490
            read 170, ians
            if (ians .ne. 'y' ) go to 1580
c
c        Print table of back-transformed data 
c
            print 1500
            if (ip) write (10,1500) 
            print 1510
            if (ip) write (10,1510) 
            do 1530 i=1,m 
               print 1520, i, bt(i,1), by(i), bketa(i)
               if (ip) write (10,1520) i, bt(i,1), by(i), bketa(i)
1530        continue
c
      end if
c
1490  format (/' Do you wish to print table of ',
     1                              'back-transformed values ?  '/)
1500  format ('1'/' observed and predicted values back-transformed'/) 
1510  format (' ','index',2x,'ind var 1',4x,'obs val',6x,'pred val'/) 
1520  format (' ',i3,1p3e13.5)
c
c------------------------- Run the statistical package -------------------------
c
1580  continue
c
      nkg2 = n + kg + 2
c
      mo2   = m / 2
      mo2p1 = mo2 + 1
c
c    -------------------------------------------------------------------
      call STAT2 ( m, n, nkg2, kg, y, weight, c, alf, a, 
c
     *            iv, t, eta, r, 
c
     *            tlabel, tlabel, ylabel, ylabel,
c
     *            clab, alflab, corlab )
c    -------------------------------------------------------------------
c  
c--  Print the correlation matrix
c
      do 1260 i=1,npkg
         do 1260 j=i,npkg
            cr(i,j) = a(i,j) / ( sqrt( a(i,i) ) * sqrt( a(j,j) ) )
            cr(j,i) = cr(i,j)
 1260 continue
c-----      print 1270
c-----1270  format ('1'/' '/' '/15x,' Correlation matrix')
c-----      if (ip) write (10,1270) 
c
c-----      call ARRAYD ( cr, npkg, npkg, maxn+maxkg )
c
c    --------------------------------------------------
      call PRCORR ( cr, maxn+maxkg, npkg, corlab, ip )
c    --------------------------------------------------  This subroutine is
c                                                           in the file STAT2.f      
c
      if ( .not. ts ) then
                      if (ip) write (10,1800)
c
            go to 770
c
      end if
c
c------------------- Fourier spectral analysis of residuals --------------------
c
1900  continue
c
      do 2000 iperd = 1, 2
c                         --- Beginning to loop to do periodogram analyses on 
c                                  the two sets of residual separately
c
      if ( iperd .eq. 1 ) then
                               ip1 = mo2 + 1
                               wavlab = 'Cosine Wave '
               else
                    ip1 = 1
                    wavlab = ' Sine Wave  '
      end if
c
      mo2 = m / 2
c
      delt =  ( t(mo2,1) - t(1,1) ) / float(mo2 - 1)
c
c   ------------------------------------------------------------------
      call PERDGR ( r(ip1), mo2, delt, nf, spwork, freq, perdg, amp,
c
     1              sump, pave, pmax, fpmax, fisher, pmsig,
c
     2              cumper, nout, flb, clb, fub, cub, whn, info )
c   ------------------------------------------------------------------
c
      if ( info .gt. 0 ) then
c
             print *, ' Error in call to PERDGR,  info =', info,
     1                '     iperd =', iperd
             if (ip) write (10, 1980) info, iperd
c
             go to 770
c
      end if
c
 1980 format('1'/' Error in call to PERDGR,  info =', i3,
     1           '     iperd =', i3 )
c
      xpl(1) = flb(1)
      ypl(1) = clb(1)
      xpl(2) = flb(2)
      ypl(2) = clb(2)
      xpl(3) = flb(2)
      ypl(3) = 1.0
      xpl(4) = fub(2)
      ypl(4) = cub(2)
      xpl(5) = fub(1)
      ypl(5) = cub(1)
c
      permx = 1.0 / fpmax
      ratio = pmax / pave
c
      if (ip) write (10,1935) wavlab, sump, pave
                  print 1935, wavlab, sump, pave
 1935 format ('1'/' FOR ', a12, 'PERIODOGRAM'/' '/
     1            ' Sum of periodogram ordinates =', 1pe12.4 /' '/
     1            ' Average periodogram ordinate =', 1pe12.4 /' '/
     1            ' Maximum periodogram ordinate ')
      if (ip) write (10,1945) fpmax, permx, pmax
                  print 1945, fpmax, permx, pmax
 1945 format (    '      freq. , period , power  = ', 
     1                    f10.6, f12.4, 1pe12.4 )
      if (ip) write (10,1940) ratio, fisher
                  print 1940, ratio, fisher
 1940 format (' '/' P_maximum / P_average  =', 1pe12.4 /
     1        ' '/'       Fisher statistic =', 1pe12.4 )
      if ( pmsig ) then 
                        write (10, 1942)
                        print 1942
              else
                   write (10,1943)
                   print 1943
      end if
 1942 format(' '/' The maximum peak is significantly above ',
     1                                      'the average value' / )
 1943 format(' '/' The maximum peak is not significantly ',
     1                                'above the average value' / )
c
      if ( whn .eq. .true. ) then
                                 if (ip) write (10, 1985) nout, nf
                                              print 1985, nout, nf
                 else 
                     write (10, 1986) nout, nf
                          print 1986, nout, nf        
      end if
 1985 format(' '/' Time-series not significantly different',
     1                  ' from white noise'/' '/' ', i6, ' / ', i6,
     2                  '  periodogram ordinates outside white',
     3                  ' noise band' )
 1986 format(' '/' Time-series differs significantly',
     1                  ' from white noise'/' '/' ', i6, ' / ', i6,
     2                  '  periodogram ordinates outside white',
     3                  ' noise band' )
c
      do 1990 i=1,nf
         spwork(i) = 0.0
 1990 continue
c
c    -----------------------------------------------------------------
      call SCPLT1 ( freq, perdg, nf, 0, 1,
c
     1              wavlab//'Residual Periodogram', 
c
     2              'Frequency', 'Variance' )
c   ------------------------------------------------------------------
      call PLOT3 ( freq, cumper, nf, 0, 1, 
c
     1             flb,       clb,        2,  0, 1,
c
     2             fub,       cub,        2,  0, 1,
c
     3             wavlab//'Cumulative Periodogram', 
c
     4             'Frequency', 'Cumulative Fraction' )
c   -----------------------------------------------------------------
      call TWOPLT ( freq, perdg, nf, freq, spwork, nf, 0, 1, 
     1              wavlab//'Residual Periodogram', 
     1              'Frequency', 'Variance', 
     2              freq, cumper, nf, xpl, ypl, 5, 0, 1,
     3              wavlab//'Cumulative Periodogram', 
     4              'Frequency', 'Cumulative Fraction', ts )
c   -----------------------------------------------------------------
      call TWOPLT ( freq, perdg, nf, freq, spwork, nf, 0, 1, 
     1              wavlab//'Residual Periodogram', 
     1              'Frequency', 'Variance', 
     2              freq, amp, nf, freq, spwork, nf, 0, 1,
     3              wavlab//'Residual Amplitude Spectrum', 
     4              'Frequency', ylabel, ts )
c   ----------------------------------------------------------------
c
c
 2000 continue
c
c          --- End of loop to do periodogram analyses on two halves separately
c
c-- Calculate periodogram of complex data obtained by combining the two waves --
c
      mo2 = m / 2
c
      do i = 1, mo2
c----         cmplxts(i) = CMPLX( r(i), r(mo2+i) )
         cmplxts(i) = CMPLX( r(mo2+i), r(i) )
      end do
c
      iou = 10
c
c    --------------------------------------------------------------------
      call CMPLXPER ( t, cmplxts, mo2, delt, nf, freq, perdg, amp,
c
     1                spwork, spwork(nf+1), ylabel, ts,
c
     1                cwrk1, cwrk1(nf+1), cwrk1(2*nf+1), cwrk1(3*nf+1),         
c
     2                cwrk2, cwrk2(mmax+1), cwrk2(2*mmax+1),
c
     3                      cwrk2(3*mmax+1), cwrk2(4*mmax+1), iou,
c
     4                       'Residual ' )
c    --------------------------------------------------------------------
c
      go to 770
c
c-------------------------------------------------------------------------------
c
1590  continue
c
c    -------------------------------------------------------
      call CH ( icode, mmax, m, kg, iv, t, y, weight, alf ) 
c    -------------------------------------------------------
c                                       subroutine to change  
      go to 770
c
c
1600  continue
c
c    --------------------------------------------------------------
      call PR ( 99999, icode, mmax, m, n, kg, p, iv, t, y, weight,
c
     *                                                       alf ) 
c    --------------------------------------------------------------
c                                               subroutine to print 
      go to 770
c
c
1610  continue
c
      print 1620
1620  format (' command code was invalid'//)
      go to 770
c
c
c------------- Save residuals and continue ------------------------
c
1630  continue
      print 1640
1640  format (' type name of file for storing residual data'/
     1        ' (maximum of 26 characters):  ')
1641  read 660, fname
      inquire (file=fname,exist=ex)
      if (.not.ex) goto 1643
      print 1642,fname
1642  format (' '/' ',a26,' already exists.  use another name.')
      goto 1641
1643  open(unit=3,file=fname,status='new')
      do 1660 i=1,m 
         write (3,1650) (t(i,j),j=1,iv),r(i)
1650     format (7(e14.8,2x)) 
1660  continue
      rewind 3
      close(unit=3,status='keep')
      if (ip) write (10,1670) fname
1670  format (//' residual data is stored in file:  ',a26)
      go to 770
c
c--------------- save current values of data ---------------------
c
1680  continue
      print 1690
1690  format (' type name of file for storing data'/
     2        ' (maximum of 26 characters):   ')
1691  read 660, fname
      inquire (file=fname, exist=ex)
      if (.not.ex) goto 1693
      print 1692,fname
1692  format (' '/' ',a26,' already exists.  use another name.')
      goto 1691
1693  open(unit=3,file=fname,status='new')
      write (3,1700) m,n,kg,iv
1700  format (' data pts: ',i4/' linear parameters: ',i4/
     2        ' nonlinear parameters: ',i4/' independent variables: ',
     3        i4/' indep. var., dep. var. & weights: ')
      do 1710 i=1,m 
1710     write (3,340) (t(i,j),j=1,iv),y(i),weight(i)
      write (3,1720) (alf(i),i=1,kg)
1720  format (' nonlinear parameters: '/ 10e14.8 )
      do 1760 i=1,kg
         do 1730 k=1,np1
1730        ia(k)=0.
         k=1
         do 1740 j=1,np1
            if (inc(i,j).ne.1) go to 1740
            ia(k)=j 
            k=k+1
1740     continue
         write (3,1750) i,(ia(k),k=1,np1)
1750     format (' Terms with nonlin parm',i3,': ',26i3,/
     1           '                       ',    5x, 26i3,/
     2           '                       ',    5x, 26i3,/
     3           '                       ',    5x, 22i3 )
1760  continue
      write (3,1770) (ndi(i),i=1,p)
1770  format (' numerical derivative increments: ' / 10e14.8)
      rewind 3
      close(unit=3,status='keep')
      if (ip) write (10,1780) fname
1780  format (//' data is stored in file:  ',a26/)
      if (.not.wrnbef) go to 790
      if (icode.eq.'sc') go to 770
c
c-------------------------- stop ---------------------------------------
c
1790  continue
      if (ip) write (10,1800) 
1800  format ('1')
c
c
      if ( ip ) then
                    rewind 10
                    close ( unit=10 , status = 'keep' )
      end if
c
      stop
c
      end 
