c
c
c*************************************************************************
c
c
      subroutine STAT2 ( m, n, nkg2, kg, y, weight, c , alf, a,
c
     *                  iv, t, eta, r, 
c
     *                  tlabel, tlabld, ylabel, ylabld,
c
     *                   clab, alflab, corlab )
c
c------------------------------------------------------------------------
c
c   provide the indicated information in response to
c   the rs command: 
c  (1)  standard deviation of the parameters
c  (2)  t-ratio
c  (3)  residual sum of squares
c  (4)  residual standard error
c  (5)  coefficient of determination
c  (6)  standardized residuals
c  (7)  plot of standardized residuals vs indendent
c       variables
c  (8)  plot of standardized residuals vs predicted
c       values
c  (9)  normal probability plot.
c
c-------------------------------------------------------------------------------
c
      PARAMETER ( mmax=3000, maxn=100, maxkg=50, miv=6 )
c
      character*1 ians
c
      character*2  kw, kl, knl
c
      character*4  corlab(maxn+maxkg)
c
      character*7  clab(maxn), alflab(maxkg)
c
      character*8  Ajlab, Bjlab, Cjlab, Djlab
c
      character*47  tlabel(miv), ylabel
      character*48  tlabld(miv), ylabld
c
      real  y(m), c(maxn), alf(maxkg), a(mmax,nkg2)
c
      real  t(mmax,iv), weight(m), eta(m), r(m), temp(mmax)
c
      real  bm(mmax,maxn+maxkg), xx1(mmax,maxn+maxkg), ph(maxn+1)
c
      real  syh(mmax), sr(mmax), srv(mmax) 
c
      logical  ip, ts, idp, idcp, itp
c
      dimension  ior(mmax)
c
c
      common /blokip/ ip
      common /blkph/ ph
      common /blkbm/ bm
      common /blkts/ ts
      common /blkitp/ itp
      common /blkidp/ idp
      common /blidcp/ idcp
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
      dimension  Asto(maxkg), Bsto(maxkg), Csto(maxkg), absC(maxkg)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c--  Initialize values
c
      npkg=n+kg
      rms=a(2,nkg2) 
      pi = acos( -1.0 )
c
c----------  Create the parameter labels  -----------
c
      if ( ncon .gt. 0 ) then
                         icon = 1
                         clab(1) = ' E( 1) '
                         clab(2) = ' F( 1) '
                         corlab(1) = 'E 01'
                         corlab(2) = 'F 01'
               else
                    icon = 0
      end if
c
      do kwav = 1, ncwav
c
         write ( kw, '(i2)' ) (icon + kwav)
c
         kc = ncon + 2*kwav - 1
         write ( kl,    '(i2)' ) (icon + kwav)
         clab(kc)   = ' E(' // kl // ') '
         clab(kc+1) = ' F(' // kl // ') '
         corlab(kc)   = 'E ' // kl
         corlab(kc+1) = 'F ' // kl
c
         kalf = 2*kwav - 1
         write ( knl,   '(i2)' ) (icon + kwav)
         alflab(kalf)   = ' B(' // knl // ') '
         alflab(kalf+1) = ' C(' // knl // ') '
         corlab(n+kalf)   = 'B ' // knl
         corlab(n+kalf+1) = 'C ' // knl
c
      end do
c
      if ( npexp .gt. 0 ) then
c
                     kc = ncon + 2*ncwav + 1
                     write ( kl, '(i2)' ) (icon + ncwav + 1)
                     clab(kc) =   ' E(' // kl // ') '
                     clab(kc+1) = ' F(' // kl // ') '
                     corlab(kc)   = 'E ' // kl
                     corlab(kc+1) = 'F ' // kl
c
                     kalf = 2*ncwav + 1
                     write ( knl, '(i2)' ) (icon + ncwav + 1)
                     alflab(kalf) = ' B(' // knl // ') '
                     corlab(n+kalf) = 'B ' // knl
c
      end if
c 
      if ( npsin .gt. 0 ) then
c
             do ksin = 1, npsin
c
                kc = ncon + 2*ncwav + 2*npexp + (2*ksin - 1)
                write ( kl, '(i2)' ) (icon + ncwav + npsin + ksin)
                clab(kc)   = ' E(' // kl // ') '
                clab(kc+1) = ' F(' // kl // ') '
                corlab(kc)   = 'E ' // kl
                corlab(kc+1) = 'F ' // kl
c
                kalf = 2*ncwav + npsin + ksin
                write( knl, '(i2)' ) (icon + ncwav + npexp + ksin)
                alflab(kalf) = ' C(' // knl // ') '
                corlab(n+kalf) = 'C ' // knl
c
             end do
c
      end if
c
c--  Print parameters, their std dev & t-ratios
c
      print 80
      if (ip) write (10,80)
80    format ('1'/' '/4x, 'Param.', 4x, 
     1            ' Estimate',6x,'Std. Deviation',5x,'t-ratio'/' ')
c
      if ( n .eq. 0 ) go to 110
c
      do 100 i=1,n
         ai=abs(a(i,i))
         bb=sqrt(ai)
         tr=c(i)/bb 
         print 90, clab(i), c(i), bb, tr 
         if (ip) write (10,90) clab(i), c(i), bb, tr
         if ( 2*(i / 2) .eq. i ) then
                                 print 105
                                 if (ip) write (10,105)
         end if
100   continue
c
   90 format (' ',4x, a7, 1p3e16.6)
  105 format( ' ' )
c
c
110   if (kg.eq.0) go to 130
c
      print 115
      if (ip) write (10,115)
  115 format (' '/' ', 4x, '----------------------------------------',
     1                     '---------------' / ' ' )
c
      do 120 i=1,kg 
         j = n + i
         aj=abs(a(j,j))
         bb=sqrt(aj)
         tr=alf(i)/bb
         print 90, alflab(i), alf(i), bb, tr
         if (ip) write (10,90) alflab(i), alf(i), bb, tr
         if ( 2*(i / 2) .eq. i ) then
                                 print 105
                                 if (ip) write (10,105)
         end if
  120 continue
c
c--  Print residual sum of squares, res mean square, res std error, and
c          coefficient of determination.
c
 130  continue
c
      rss=(a(1,nkg2))**2
      rsd=sqrt(abs(rms))
      print 140, rss
      if (ip) write (10,140) rss
140   format (' '/' '/10x, ' wtd residual sum of squares:  ',1pe16.6)
      print 150, rms
      if (ip) write (10,150) rms
150   format (' '/' ',10x,'wtd residual mean square:  ',3x,1pe16.6)
      print 160, rsd
      if (ip) write (10,160) rsd
160   format (' '/' ',10x,'wtd residual standard error:  ',1pe16.6/' ')
      yb=0.
      sws=0.
      do 170 i=1,m
         sws=sws+weight(i)
170      yb=yb+y(i)*weight(i) 
      yb=yb/sws
      ctss=0.
      do 180 i=1,m
180      ctss=ctss+weight(i)*(y(i)-yb)**2
      r2=1.-rss/ctss
      print 190, r2 
      if (ip) write (10,190) r2
190   format (' coefficient of determination (r-square): ',1pe16.6/' ')
c
c
c--  Convert coefficients of sine and cosine terms (E and F) to amplitudes
c         and phases (A and D) -- and also their uncertainties
c
      write (10,192)
  192 format( '1')
c
      write (10,193)
  193 format (' '/' '/' For Damped Sine/Cosine Pairs '/
     1                ' ---------------------------- ' )
c
      nlp = ncon + 2*ncwav + 2*npexp + 2*npsin
c
      do kwav = 1, ncwav
c
         write ( kw, '(i2)' ) (icon + kwav)
c
         Ajlab = ' A(' // kw // ') ='
         Bjlab = ' B(' // kw // ') ='
         Cjlab = ' C(' // kw // ') ='
         Djlab = ' D(' // kw // ') ='
c
         kc   = ncon + 2*kwav - 1
         kalf =        2*kwav - 1
c
         Ec    = c(kc)
         sigE  = SQRT( a(kc,kc) )
         Fc    = c(kc+1)
         sigF  = SQRT( a(kc+1,kc+1) )
         covEF = a(kc,kc+1)
c
         Bc = alf(kalf)
         Cc = alf(kalf+1)
c
         kalf = nlp + kalf
c
         if ( kg .gt. 0 ) then
                              sigB = SQRT( a(kalf,kalf) )
                              sigC = SQRT( a(kalf+1,kalf+1) )
                  else
                       sigB = 0.0
                       sigC = 0.0
         end if
c
         tB = Bc / sigB
         tC = Cc / sigC
c
c   ----------------------------------------------------------------
         call EFtoAB ( Cc, Ec, sigE, Fc, sigF, covEF, 
c
     1                                         Ac, sigA, Dc, sigD )
c   ----------------------------------------------------------------
c
         tA = Ac / sigA
c
         write (10,195) Ajlab, Ac, sigA, tA,
     1                  Bjlab, Bc, sigB, tB,
     2                  Cjlab, Cc, sigC, tc,
     3                  Djlab, Dc, sigD
c
         Asto(kwav) = Ac
         Bsto(kwav) = Bc
         Csto(kwav) = Cc
         absC(kwav) = abs( Cc )
c
      end do
c
  195 format (' ',
     *  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [V],     t =', 1pe9.2,
     1  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [ms-1],  t =', 1pe9.2,
     2  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [kHz],   t =', 1pe9.2,
     3  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [ms]' )
c
c
      if ( npexp .eq. 1 ) then
c
         write (10,194)
  194    format(' '/' '/' For Zero-frequency Sine/Cosine Pair ',
     1                  '(Pure Exponential Decay) '
     2                 /' ------------------------------------',
     3                  '------------------------ ' )
c
         write ( kw, '(i2)' ) (icon + ncwav + 1)
c
         Ajlab = ' A(' // kw // ') ='
         Bjlab = ' B(' // kw // ') ='
         Djlab = ' D(' // kw // ') ='
c
         kc   = ncon + 2*ncwav + 1
         kalf =        2*ncwav + 1
c
         Ec    = c(kc)
         sigE  = SQRT( a(kc,kc) )
         Fc    = c(kc+1)
         sigF  = SQRT( a(kc+1,kc+1) )
         covEF = a(kc,kc+1)
c
         Bc = alf(kalf)
c
         kalf = nlp + kalf
c
         if ( kg .gt. 0 ) then
                              sigB = SQRT( a(kalf,kalf) )
                   else
                        sigB = 0.0
         end if
c
         tB = Bc / sigB
c
         Cc = 0.0
c   ----------------------------------------------------------------
         call EFtoAB ( Cc, Ec, sigE, Fc, sigF, covEF, 
c
     1                                         Ac, sigA, Dc, sigD )
c   ----------------------------------------------------------------
c
         tA = Ac / sigA
c
         write (10,196) Ajlab, Ac, sigA, tA,
     1                  Bjlab, Bc, sigB, tB,
     3                  Djlab, Dc, sigD

c
      end if
c
  196 format (' ',
     *  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [V],     t =', 1pe9.2,
     1  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [ms-1],  t =', 1pe9.2,
     3  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [rad]' )
c
c
      if ( npsin .ne. 0 ) then
c
         write (10,197)
  197    format(' '/' '/' For Undamped Sine/Cosine Pairs '/
     1                  ' ------------------------------ ' )
c
         do kwav = 1, npsin
c
            write ( kw, '(i2)' ) (icon + ncwav + npexp + kwav)
c
            Ajlab = ' A(' // kw // ') ='
            Cjlab = ' C(' // kw // ') ='
            Djlab = ' D(' // kw // ') ='
c
            kc   = ncon + 2*ncwav + 2*npexp + (2*kwav - 1)
            kalf =        2*ncwav +   npexp +    kwav
c
            Ec    = c(kc)
            sigE  = SQRT( a(kc,kc) )
            Fc    = c(kc+1)
            sigF  = SQRT( a(kc+1,kc+1) )
            covEF = a(kc,kc+1)
c
            Cc = alf(kalf)
c
            kalf = nlp + kalf
c
            if ( kg .gt. 0 ) then
                                 sigC = SQRT( a(kalf,kalf) )
                      else
                           sigC = 0.0
            end if
c
            tC = Cc / sigC
c
c   ----------------------------------------------------------------
            call EFtoAB ( Cc, Ec, sigE, Fc, sigF, covEF, 
c
     1                                         Ac, sigA, Dc, sigD )
c   ----------------------------------------------------------------
c
            tA = Ac / sigA
c
            write (10,198) Ajlab, Ac, sigA, tA,
     2                     Cjlab, Cc, sigC, tc,
     3                     Djlab, Dc, sigD

c
         end do
c
      end if
c
  198 format (' ',
     *  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [V],     t =', 1pe9.2,
     2  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [kHz],   t =', 1pe9.2,
     3  /' ', a8, 1pe15.7, ' +/-', 1pe15.7, ' [ms]' )
c
c
c-- Check for correlations between amplitudes, frequencies and damping rates
c
      if ( ncwav .ge. 3 ) then
c
c    -----------------------------------------------------------------
             call SCPLT1 ( Csto, Asto, ncwav, 1, 0, 
c
     1                     'Amplitude vs. Frequency',  
c
     2                                          'C [kHz]', 'A [V]' )
c    -----------------------------------------------------------------
             call SCPLT1 ( absC, Asto, ncwav, 1, 0, 
c
     1                     'Amplitude vs. ||Frequency||',  
c
     2                                      '||C|| [kHz]', 'A [V]' )
c    -----------------------------------------------------------------
             call SCPLT1 ( Csto, Bsto, ncwav, 1, 0, 
c
     1                     'Damping Constant vs. Frequency',  
c
     2                                      'C [kHz]', 'B [ms^-^1]' )                                      
c    -----------------------------------------------------------------
             call SCPLT1 ( absC, Bsto, ncwav, 1, 0, 
c
     1                     'Damping Constant vs. ||Frequency||',  
c
     2                                  '||C|| [kHz]', 'B [ms^-^1]' )                                      
c    -----------------------------------------------------------------
             call SCPLT1 ( Bsto, Asto, ncwav, 1, 0,
c
     1                     'Amplitude vs. Damping Constant',
c
     2                                        'B [ms^-^1]', 'A [V]' )
c    ------------------------------------------------------------------
c      
      end if
c
c
c-------------------------------------------------------------------------------
c
c--  Calculate standardized residuals
c
200   if (n.eq.0) go to 220
      do 210 i=1,m
         ii=i
c
c       -----------------------------------------------------
         call CALC ( ii, kg, n, mmax, iv, c, alf, t, dummy )
c       -----------------------------------------------------
c
         do 210 j=1,n
210         bm(i,j)=ph(j)
220   do 240 j=1,npkg
         do 240 i=1,m
            sum=0.
            do 230 k=1,npkg
230            sum=sum+bm(i,k)*a(k,j)
240         xx1(i,j)=sum
      do 260 i=1,m
         sum=0.
         do 250 k=1,npkg
250         sum=sum+xx1(i,k)*bm(i,k)
         syh(i)=sqrt(abs(sum))
260      sr(i)=r(i)/sqrt(abs(rms/weight(i)-sum))
c
c-- Calculate expected standardized residuals ( for normal probability plot)
c
      if (m.gt.10) go to 410
c
      do 400 i=1,m
         pp = (3.0*i - 1.0) / (3.0*m + 1.0)
         syh(i) = 4.91*( pp**(0.14) - (1.0 - pp)**(0.14) )
400   continue
c
      go to 430
c
410   continue
c      
      pp = 2.0 / (3.0*m + 1.0 )
      syh(1) = 4.91*( pp**(0.14) - (1.0 - pp)**(0.14) )
c
      pi8 = pi / 8
      pi4 = pi / 4
c
      do 420 i=2,m
         pp = (i - pi8) / ( m + 1.0 - pi4 )
         syh(i) = 4.91*( pp**(0.14) - (1.0 - pp)**(0.14) )
420   continue
c
      mo2 = m / 2
      mo2p1 = mo2 + 1
c
c    -----------------------------
430   call ORPLOT (sr, m, ior )
c    -----------------------------
      do 440 i=1,m
         srv(i)=sr(ior(i))
  440 continue
c
c
c--  VG plot of (1) std. residuals vs. ind. var. 
c               (2) std. residuals vs. dep. var.
c               (3) obs. std. resids. vs. expected std. resids. 
c
      do 360 i=1,iv 
c
         do 355 j=1,m
            temp(j) = 0.0
  355    continue
c
      mo2 = m / 2
      mo2p1 = mo2 + 1
c
c         if (ts.and.itp) then
c    --------------------------------------------------------------------
c                   call SCPLT2 ( t(1,i), sr, mo2, 
c
c     *                           t(1,i), temp, mo2, 0, 1,
c
c     *                           'Sine Wave Std Resids vs Indep Var',
c
c     *                                       tlabel(i), 'Std Resids' )
c------------------------------------------------------------------------
c                   call SCPLT2 ( t(mo2p1,i), sr(mo2p1), mo2, 
c
c     *                           t(mo2p1,i), temp(mo2p1), mo2, 0, 1,
c
c     *                           'Cosine Wave Std Resids vs Indep Var',
c
c     *                                        tlabel(i), 'Std Resids' )
c------------------------------------------------------------------------
c         end if
c
c         if (.not.ts) then
c    --------------------------------------------------------------------
c               call SCPLT2 ( t(1,i), sr, mo2, 
c
c     *                       t(1,i), temp, mo2, 1, 1,
c
c     *                       'Sine Wave Std Resids vs Indep Var',
c
c     *                                    tlabel(i), 'Std Resids' )
c    --------------------------------------------------------------------
c               call SCPLT2 ( t(mo2p1,i), sr(mo2p1), mo2, 
c
c     *                       t(mo2p1,i), temp(mo2p1), mo2, 1, 1,
c
c     *                       'Cosine Wave Std Resids vs Indep Var',
c
c     *                                    tlabel(i), 'Std Resids' )
c    --------------------------------------------------------------------
c      end if
c
c
         if (ts.and.itp) then
c   ---------------------------------------------------------------------
                call THRPLT ( t(1,i), sr, mo2, 
c
     *                        t(1,i), temp, mo2, 0, 1, 
c
     *                        'Sine Wave Std. Resids. vs. Ind. Var.',
c
     *                        tlabel(i), 'Std. Resids.',
c
     *                        eta, sr, mo2, 
c
     *                        eta, temp, mo2, 17, 1,
c
     *                        'Std Resids vs Dep Var',
c
     *                        ylabel, 'Std. Resids.',
c
     *                        syh, srv, m, 17, 0, 
c
     *                        'Normal Prob. Plot',
c
     *              'Expected Std. Resids.', 'Obs. Std. Resids.', ts )
c   ----------------------------------------------------------------------
                call THRPLT ( t(mo2p1,i), sr(mo2p1), mo2, 
c
     *                        t(mo2p1,i), temp(mo2p1), mo2, 0, 1, 
c
     *                        'Cosine Wave Std. Resids. vs. Ind. Var.',
c
     *                        tlabel(i), 'Std. Resids.',
c
     *                        eta(mo2p1), sr(mo2p1), mo2, 
c
     *                        eta(mo2p1), temp(mo2p1), mo2, 17, 1,
c
     *                        'Std Resids vs Dep Var',
c
     *                        ylabel, 'Std. Resids.',
c
     *                        syh, srv, m, 17, 0, 
c
     *                        'Normal Prob. Plot',
c
     *              'Expected Std. Resids.', 'Obs. Std. Resids.', ts )
c   ----------------------------------------------------------------------
         end if
c
         if (.not.ts) then
c   ----------------------------------------------------------------------      
               call THRPLT ( t(1,i), sr, mo2, 
c
     *                       t(1,i), temp, mo2, 17, 1,
c
     *                       'Sine Wave Std. Resids. vs. Ind. Var.',
c
     *                       tlabel(i), 'Std. Resids.',
c
     *                       eta, sr, mo2, 
c
     *                       eta, temp, mo2, 17, 1,
c
     *                       'Std Resids vs Dep Var',
c
     *                       ylabel, 'Std. Resids.',
c
     *                       syh, srv, m, 17, 0,
c
     *                   'Normal Prob. Plot',
c
     *              'Expected Std. Resids.', 'Obs. Std. Resids.', ts )         
c   ---------------------------------------------------------------------
               call THRPLT ( t(mo2p1,i), sr(mo2p1), mo2, 
c
     *                       t(mo2p1,i), temp(mo2p1), mo2, 17, 1,
c
     *                       'Cosine Wave Std. Resids. vs. Ind. Var.',
c
     *                       tlabel(i), 'Std. Resids.',
c
     *                       eta(mo2p1), sr(mo2p1), mo2, 
c
     *                       eta(mo2p1), temp(mo2p1), mo2, 17, 1,
c
     *                       'Std Resids vs Dep Var',
c
     *                       ylabel, 'Std. Resids.',
c
     *                       syh, srv, m, 17, 0,
c
     *                   'Normal Prob. Plot',
c
     *              'Expected Std. Resids.', 'Obs. Std. Resids.', ts )         
c   ---------------------------------------------------------------------
         end if
c
c    -----------------------------------------------------------
c         call SCPLT2 ( eta, sr, mo2, 
c     *                 eta, temp, mo2, 1, 1,
c     *                'Sine Wave Std. Resids. vs. Dep. Var.',
c     *                 ylabel, 'Std. Resids.' )
c    ------------------------------------------------------------
c         call SCPLT2 ( eta(mo2p1), sr(mo2p1), mo2, 
c     *                 eta(mo2p1), temp(mo2p1), mo2, 1, 1,
c     *                 'Cosine Wave Std. Resids. vs. Dep. Var.',
c     *                 ylabel, 'Std. Resids.' )
c    ------------------------------------------------------------
c
c    ------------------------------------------------------------------
c         call SCPLT1 ( syh, srv, m, 1, 0, 
c     *                 'Normal Probability Plot',
c     *                 'Expected Std. Resids.', 'Obs. Std. Resids.' )
c    ------------------------------------------------------------------
c
c
360   continue
c
c
c-------------------------------------------------------------------------------
c
c--  Print table of obs. & pred. values along with std. res. if desired
c
      print 40
40    format (' Print table of obs, pred val, resid & std. resid?  ')
      read (5,20) ians
20    format ( a1 )
      if ( (ians .ne. 'y') .and. (ians .ne. 'Y') ) go to 330 
c
      do 300 inv=1,iv
         print 270, inv
         if (ip) write (10,270) inv
         do 290 i=1,m
            print 280, i,t(i,inv),y(i),eta(i),r(i),sr(i)
            if (ip) write (10,280) i,t(i,inv),y(i),eta(i),r(i),sr(i)
290      continue 
300   continue
270   format ('1'/' index',2x,'ind var', i2, 4x,'obs val',5x,'pred val',
     2             6x, 'resid',8x,'std res'/)
280   format (' ',i3,1p5e13.5)
c
c-------------------------------------------------------------------------------
 330  continue
c
c    ++++++++
 460  return
c    ++++++++
c
      end 
c
c
c*******************************************************************************
c
c
      subroutine PRCORR ( cr, ldcr, npar, corlab, ip )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c     Subroutine for printing the correlation matrix cr
c
      dimension  cr(ldcr,npar)
c
      character*4  corlab(npar)
c
      logical  ip
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      do j = 1, npar, 16
c
         minj = j
         maxj = min( (j + 15) , npar )
                    print 20, ( corlab(jj), jj = minj, maxj )
         if (ip) write(10,20) ( corlab(jj), jj = minj, maxj )
c
         do i = 1, npar
c
            if ( (i .eq. 41) .or. (i .eq. 81) .or. (i .eq. 121) ) then
                              print 20, ( corlab(jj), jj = minj, maxj )
                   if (ip) write(10,20) ( corlab(jj), jj = minj, maxj )
             end if
c
                      print 40, corlab(i), ( cr(i,jj), jj = minj, maxj )
           if (ip) write(10,40) corlab(i), ( cr(i,jj), jj = minj, maxj )
c         
         end do
c
      end do
c
   20 format('1'/' '/' Correlation Matrix',
     1          /' '/' ',5x, 16( 3x, a4 ) / )
c
   40 format(' ', a4, ' ', 0p16f7.3 )
c
c    ++++++++
      return
c    ++++++++
c
      end      
      
