      subroutine errstat(modvis2,chisqm)
c
c     Calculate and display varios statistics on data and model
c     It produces six plots.
c
c     C. A. Hummel      30 October 1992
c
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'
c
      real*8 modvis2(mxrec,mxfilt)
      real*8 mhours(mxpl),mh,lambda,dlambda,mv,mp,u,v,seprad,parad
      real*4 y(mxscn),t(mxscn)
      real*4 timmax,timmin,timstep
      integer*4 firstid,lastid
      logical query,hardcopy,doed
      character header*80,device*4,plotfile*16,command*1
      
      integer*4 ma
      parameter (ma=7)
      integer*4 lista(ma)
      real*8 p(ma),pa(ma),covar(ma,ma)
      real*8 y8(mxscn),t8(mxscn),s8(mxscn),chisq,chisqm
      real*4 x4data(mxrec2),y4data(mxrec2),e4data(mxrec2),dev2(mxrec2)
      real*4 calerr(mxrec,mxfilt)
      real*4 dev(mxrec)
      real*4 x1(40),y1(40),x2(40),y2(40),autocorr(12),npairs(12),lag(12)
 
      external fpoly
c
c     First, do statistics of observed-model data
c
    1 write(outc,*)'Plot 1: Error distribution'
      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
    2 format(a4)
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='ploterr1.lpm'//'/lpm'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='ploterr1.ps'//'/vps'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='ploterr1.ps'//'/ps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
c
      doed=.false.
      if(query(' Edit data based on normalized deviation ? (y/n) '))then
       write(outc,'(a,$)')' Enter maximum norm. deviation: '
       read(inc,*)devmax
       doed=.true.
      endif
c
c.....Do four plots of error budget
      i2bad=0
      nrec=ndata*nfilt
      chif=real(chisqm)/nfilt/ndata
      if(calib)chif=1.0
      do ig=1,nfilt
       ibad=0
       xhigh=0.
       do id=1,ndata
        if(viserr(id,ig).ne.2.D0)then
         x4data((ig-1)*ndata+id-i2bad)=real(modvis2(id,ig))
         y4data((ig-1)*ndata+id-i2bad)=real(obsvis2(id,ig))
         dev(id-ibad)=real((obsvis2(id,ig)-modvis2(id,ig))
     &                     /viserr(id,ig))
         dev2((ig-1)*ndata+id-i2bad)=dev(id-ibad)
         e4data((ig-1)*ndata+id-i2bad)=real(viserr(id,ig))
         xhigh=max(xhigh,abs(dev(id-ibad)))
         if(doed)then
          if(abs(dev(id-ibad)).gt.devmax)then
           errf=abs(dev(id-ibad))/devmax
           viserr(id,ig)=viserr(id,ig)*dble(errf)
          endif
         endif
c        if(x4data((ig-1)*ndata+id-i2bad).lt.0.0312d0.and.
c    &      x4data((ig-1)*ndata+id-i2bad).gt.0.0144d0.and.
c    &      dev2((ig-1)*ndata+id-i2bad).lt.3.2d0.and.
c    &      dev2((ig-1)*ndata+id-i2bad).gt.0.0d0)then
c         write(15,*)date(id),baseline(id),lambda0(id,ig)
c         ibad=ibad+1
c         i2bad=i2bad+1
c        endif
        else
         ibad=ibad+1
         i2bad=i2bad+1
        endif
       enddo
       xhigh=max(4.,xhigh)
       xhigh=min(8.,xhigh)
       xlow=-xhigh
       nbin=int(5.*xhigh)
       if(ig.eq.1)then
        if(device.eq.'v'.or.device.eq.'V')then
         call pgbegin(0,plotfile,2,4)
        else
         call pgbegin(0,plotfile,2,2)
        endif
        call pgsch(2.0)
        hheight=real(nrec/5)/3.
        call pgsci(1)
        call pgenv(xlow,xhigh,0.0,hheight,0,0)
        call pglabel('Normalized Deviation',' ',
     &               'Error distribution (red channel)')
        call pgsci(2)
        n=ndata-ibad
        call pghist(n,dev,xlow,xhigh,nbin,1)
        do i=1,40  
         x1(i)=xlow+(i-1)*(xhigh-xlow)/40
         if(x1(i).gt.7)then
          y1(i)=0.
          y2(i)=0.
         else
          y1(i)=1./sqrt(2.0*3.1415*chif)*exp(-0.5/chif*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
          y2(i)=1./sqrt(2.0*3.1415)*exp(-0.5*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
         endif
        enddo
        call pgsls(4)
        call pgsci(10)
        call pgline(40,x1,y2)
        call pgsci(2)
        call pgsls(1)
        call pgline(40,x1,y1)
       elseif(ig.eq.2)then
        call pgsci(1)
        call pgenv(xlow,xhigh,0.0,hheight,0,0)
        call pglabel('Normalized Deviation',' ',
     &               'Error distribution (green channel)')
        call pgsci(3)
        n=ndata-ibad
        call pghist(n,dev,xlow,xhigh,nbin,1)
        do i=1,40  
         x1(i)=xlow+(i-1)*(xhigh-xlow)/40
         if(x1(i).gt.7)then
          y1(i)=0.
          y2(i)=0.
         else
          y1(i)=1./sqrt(2.0*3.1415*chif)*exp(-0.5/chif*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
          y2(i)=1./sqrt(2.0*3.1415)*exp(-0.5*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
         endif
        enddo
        call pgsls(4)
        call pgsci(10)
        call pgline(40,x1,y2)
        call pgsci(3)
        call pgsls(1)
        call pgline(40,x1,y1)
       elseif(ig.eq.3)then
        call pgsci(1)
        call pgenv(xlow,xhigh,0.0,hheight,0,0)
        call pglabel('Normalized Deviation',' ',
     &               'Error distribution (blue channel)')
        call pgsci(11)
        n=ndata-ibad
        call pghist(n,dev,xlow,xhigh,nbin,1)
        do i=1,40  
         x1(i)=xlow+(i-1)*(xhigh-xlow)/40
         if(x1(i).gt.7)then
          y1(i)=0.
          y2(i)=0.
         else
          y1(i)=1./sqrt(2.0*3.1415*chif)*exp(-0.5/chif*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
          y2(i)=1./sqrt(2.0*3.1415)*exp(-0.5*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
         endif
        enddo
        call pgsls(4)
        call pgsci(10)
        call pgline(40,x1,y2)
        call pgsci(11)
        call pgsls(1)
        call pgline(40,x1,y1)
        n=nrec-i2bad
        call pgsci(1)
        hheight=hheight*3.
        call pgenv(xlow,xhigh,0.0,hheight,0,0)
        call pglabel('Normalized Deviation',' ',
     &               'Error distribution (all channels)')
        call pghist(n,dev2,xlow,xhigh,nbin,1)
        do i=1,40  
         x1(i)=xlow+(i-1)*(xhigh-xlow)/40
         if(x1(i).gt.7)then
          y1(i)=0.
          y2(i)=0.
         else
          y1(i)=1./sqrt(2.0*3.1415*chif)*exp(-0.5/chif*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
          y2(i)=1./sqrt(2.0*3.1415)*exp(-0.5*x1(i)**2)
     &         *real(n)/(real(nbin)/(xhigh-xlow))
         endif
        enddo
        call pgsls(4)
        call pgsci(10)
        call pgline(40,x1,y2)
        call pgsci(1)
        call pgsls(1)
        call pgline(40,x1,y1)
       endif
      enddo
      call pgend
c
c.....Now some more useful plots for binaries
      if(.not.calib)then
       write(outc,*)'Plot 2: Errors vs model visibility'
       write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
       read(inc,2)device
       nx=1
       ny=1
       if(device.eq.'h'.or.device.eq.'H')then
        plotfile='ploterr2.lpm'//'/lpm'
       else if(device.eq.'v'.or.device.eq.'V')then
        plotfile='ploterr2.ps'//'/vps'
        ny=2
       else if(device.eq.'p'.or.device.eq.'P')then
        plotfile='ploterr2.ps'//'/ps'
       else if(device.eq.'s'.or.device.eq.'S')then
        plotfile='/xs'
       else
        return
       endif
       call pgbegin(0,plotfile,nx,ny)
       call pgsch(1.5)
c .... Plot # 1
       ymax=0.
       do i=1,n
        ymax=max(ymax,abs(dev2(i)))
        x4data(i)=sqrt(x4data(i))
       enddo
       ymax=min(10.,ymax)
       ymin=-ymax
       call pgenv(0.0,1.0,ymin,ymax,0,0)
       call pglabel('Model Visibility','Normalized Deviation',' ')
       call pgpoint(n,x4data,dev2,1)
       call pgsci(3)
       x1(1)=0.
       x1(2)=1.
       y1(1)=0.
       y1(2)=0.
       call pgline(2,x1,y1)
       call pgsci(1)
c .... Plot # 2
       ymax=0
       ymin=0
       do i=1,n
        e4data(i)=e4data(i)/2./x4data(i)
        ymax=max(ymax,e4data(i))
        ymin=min(ymin,e4data(i))
       enddo
       ymax=min(0.1,ymax)
       call pgenv(0.0,1.0,ymin,ymax,0,0)
       call pglabel('Model Visibility','Visibility uncertainty',' ')
       call pgpoint(n,x4data,e4data,1)
c .... Plot # 3
c      init=0
c      do ig=1,nfilt
c       do id=1,ndata
c        y4data(id)=real(viserr(id,ig))
c       enddo
c       if(ig.eq.1)call pgsci(2)
c       if(ig.eq.2)call pgsci(3)
c       if(ig.eq.3)call pgsci(11)
c       call pghist(ndata,y4data,ymin,ymax,20,init)
c       init=1
c      enddo
c      call pgsci(1)
c      call pglabel('Vis. formal error','N','Three channels')
       call pgend
      endif
c
c     Now do analysis on night by night basis
c
      write(outc,*)'Plot 3: Errors on nightly basis'
      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='ploterr3.lpm'//'/lpm'
       iy=4
       csize=2.5
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='ploterr3.ps'//'/ps'
       iy=4
       csize=2.5
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='ploterr3.ps'//'/vps'
       iy=4
       csize=2.5
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
       iy=3
       csize=2.5
      else
       return
      endif
c
      call pgbegin(0,plotfile,3,iy)
      call pgsch(csize)
      firstid=1
      lastid=1
      timmin=7
      timmax=7
      timstep=0.2
      i2bad=0
      ical=0
      do i=1,12
       autocorr(i)=0.
       npairs(i)=0.
       lag(i)=real(i)
      enddo
c
c     Go through the whole data file and ...
c 
      DO ID=1,NDATA
       timmax=max(timmax,real(hours(id)))
       timmin=min(timmin,real(hours(id)))
c
c      ... output plot, when date or baseline changes
c
       if(date(id+1).ne.date(id).or.baseline(id+1).ne.baseline(id)
     &                          .or.id+1.gt.ndata)then

       timr=timmax-timmin
       timmin=timmin-timr/10.
       timmax=timmax+timr/10.
c      loop over nfilt plots
       do ig=1,nfilt
        call pgsci(5)
        ymin=-0.2
        ymax=+0.2
        call pgenv(timmin,timmax,ymin,ymax,0,1)
        write(header,3)
     &  starno,date(id),baseline(id),int(lambda0(id,ig))
    3   format('FK#',i4,1x,i6,1x,' B.l. ',i2,1x,i3,' nm')
        if(calib)write(header(1:7),'(a)')'calibr.'
        csize=csize+1
        call pgsch(csize)
        call pglabel('','',header)
        csize=csize-1
        call pgsch(csize)
        if(ig.eq.1)call pgsci(2)
        if(ig.eq.2)call pgsci(3)
        if(ig.eq.3)call pgsci(11)
        do n=1,lastid-firstid+1
         y8(n)=obsvis2(n+firstid-1,ig)-modvis2(n+firstid-1,ig)
         s8(n)=viserr(n+firstid-1,ig)
         t8(n)=hours(n+firstid-1)
         y(n)=real(y8(n))
         t(n)=real(t8(n))
         if(s8(n).ne.2.D0)then
          call pgpoint(1,t(n),y(n),1)
          obsvis2lo=y(n)-real(s8(n))
          obsvis2hi=y(n)+real(s8(n))
          call pgerry(1,t(n),obsvis2lo,obsvis2hi,1.0)
         endif
        enddo
        call pgsci(6)
c
c ..... Do autocorrelation
        do n=1,lastid-firstid
         do 20 j=n+1,lastid-firstid+1
          do i=1,12
           if(t(j)-t(n).le.lag(i))then
            if(s8(n).ne.2.D0)then
             autocorr(i)=autocorr(i)+y(n)*y(j)
             npairs(i)=npairs(i)+1.
            endif
            goto 20
           endif
          enddo
   20    continue
        enddo
c
c ..... Do analytical fit to residuals
c       do i=1,ma
c        lista(i)=i
c        pa(i)=0.D0
c       enddo
c       k=lastid-firstid+1
c       mfit=min(7,abs(2*(k/2)-1))
c       mfit=4
c       call lfit(t8,y8,s8,k,pa,ma,lista,mfit,covar,ma,chisq,fpoly)
c
c ..... Plot analytical fit
c       do k=1,mxpl
c        mhours(k)=dble(timmin+timstep*(k-1))
c        t(k)=real(mhours(k))
c        call fpoly(mhours(k),p,mfit)
c        y(k)=0.
c        do j=1,mfit
c         y(k)=y(k)+real(p(j)*pa(j))
c        enddo
c        if(mhours(k).gt.dble(timmax))goto 40
c       enddo
c 40    ndat=min(k,mxpl)
c       call pgline(ndat,t,y)
c
c ..... Alternatively, do a running average
        do k=1,mxpl
         mhours(k)=dble(timmin+timstep*(k-1))
         t(k)=real(mhours(k))
         y(k)=0.
         swgn=0.
         do n=1,lastid-firstid+1
          if(s8(n).ne.2.D0)then
           wg=real(exp(-(t8(n)-mhours(k))**2))  ! Gaussian
           wn=1./real(s8(n)*s8(n))              ! Noise
           y(k)=y(k)+real(y8(n))*wg*wn
           swgn=swgn+wg*wn
          endif
         enddo
         if(swgn.ne.0)then
          y(k)=y(k)/swgn
         else
          y(k)=0.
         endif
         if(mhours(k).gt.dble(timmax))goto 40
        enddo
   40   ndat=min(k,mxpl)
        call pgline(ndat,t,y)
c
c ..... Fill in calibration array calerr
        icaln=0
        do k=1,ndat,4
         if(mhours(k).gt.t8(1)               .and.
     &      mhours(k).lt.t8(lastid-firstid+1))then
          if(ical.lt.mxrec)then
           ical=ical+1
           icaln=icaln+1
           m=1
           n=1
           num=1
           md=date(firstid)
           mb=baseline(firstid)
           mh=mhours(k)
           lambda=lambda0(firstid,ig)
           dlambda=dlambda0(firstid,ig)
           if(.not.calib)then
            call calcuv(n,m,1,1,starno,md,mb,mh,lambda,u,v)
            mode=0
            call binavis(n,m,1,1,md,mb,mh,lambda,dlambda,mv,mp,
     &                   u,v,seprad,parad,mode,ftcurve)
           else
            mv=1.D0
           endif
           calerr(ical,ig)=1.+y(k)/real(mv)
          else
           write(outc,*)'*** ERROR: array calerr about to overflow!'
           pause
           goto 50
          endif
         endif
        enddo 
   50   if(ig.ne.3)ical=ical-icaln
c
c ..... Fill in random error array dev2 
        do n=1,lastid-firstid+1
         if(s8(n).ne.2.D0)then
          k=(ig-1)*(lastid-firstid+1)+(n+firstid-1)-i2bad
          j=int((t8(n)-dble(timmin))/dble(timstep))+1
          if(j.gt.mxpl)j=mxpl
          dev2(k)=real((y8(n)-y(j))/s8(n))
         else
          i2bad=i2bad+1
         endif
        enddo
c
       enddo      ! End of ig-loop
c
c .... Now prepare for the next date
       timmin=7
       timmax=7
       firstid=lastid+1
       lastid=firstid-1
c
c     Continue here, if date did not change
      endif
      lastid=lastid+1
      ENDDO       ! End of ID-loop
      call pgend
c
c     Now plot autocorrelation
c
      write(outc,*)'Plot 4: Error autocorrelation'
      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='ploterr4.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='ploterr4.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='ploterr4.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
      call pgbegin(0,plotfile,1,1)
      call pgsci(1)
      call pgsch(1.5)
c
      ymin=0.
      ymax=0.
      do i=1,12
       if(npairs(i).ne.0)then
        autocorr(i)=autocorr(i)/npairs(i)
        ymin=min(ymin,autocorr(i))
        ymax=max(ymax,autocorr(i))
       else
        autocorr(i)=0.
       endif
      enddo
      call pgenv(0.,12.,ymin,ymax,0,1)
      call pglabel('Lag [h]',' ','Autocorrelation of deviation')
      call pgsci(3)
      call pgpoint(12,lag,autocorr,0)
      call pgend
c
c     Now plot histogram of the deviation of observed visibility
c     from model+polynomial fit
c
      write(outc,*)'Plot 5: True random errors'
      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='ploterr5.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='ploterr5.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='ploterr5.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
      call pgbegin(0,plotfile,1,1)
      call pgsch(1.5)
      call pgsci(1)
      n=nrec-i2bad
      call pghist(n,dev2,-4.0,4.0,20,0)
      call pglabel('Normalized Deviation from smooth residual',' ',
     &             'True Random Error Distribution')
      call pgline(40,x1,y1)        ! Theoretical curve, see above
      call pgend
c
c     Now plot histogram of the calibration error
c
      c0=real(nrec)/20.
      x0=1.0
      a0=0.1
      xmin=0.8
      xmax=1.2
c
   30 write(outc,*)'Plot 6: Calibration errors'
      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
      hardcopy=.false.
   34 if(device.eq.'h'.or.device.eq.'H')then
       plotfile='ploterr6.lpm'//'/lpm'
       hardcopy=.true.
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='ploterr6.ps'//'/ps'
       hardcopy=.true.
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='ploterr6.ps'//'/vps'
       hardcopy=.true.
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
      call pgbegin(0,plotfile,1,1)
      call pgsch(1.5)
      do ig=1,nfilt
       do i=1,ical
        dev(i)=calerr(i,ig)
        dev2((ig-1)*ical+i)=dev(i)
       enddo
       if(ig.eq.1)then
        call pgsch(1.5)
        call pgsci(1)
        nrec=ical*nfilt
        hheight=real(nrec/5)
        call pgenv(xmin,xmax,0.0,hheight,0,0)
        call pglabel('Deviation of smooth residual from model',' ',
     &               'Calibration Error Distribution')
        call pgsci(2)
       endif
       if(ig.eq.2)call pgsci(3)
       if(ig.eq.3)call pgsci(11)
       call pghist(ical,dev,xmin,xmax,20,1)
      enddo
      call pgsci(1)
      call pghist(nrec,dev2,xmin,xmax,20,1)
      do i=1,40
       x1(i)=xmin+(xmax-xmin)/40*(i-1)
       y1(i)=c0/(1.+(x1(i)-x0)**2/a0**2)
      enddo
      call pgline(40,x1,y1)
      if(hardcopy)then
       call pgend
       return
      endif
   31 write(outc,*)'Click on window and enter command:'
      call pgcurse(cx,cy,command)
      if(command.eq.'f')then
       write(outc,32)c0,a0,x0
   32  format('Enter c0, a0, x0 (ctrl z for [',
     &          1x,f4.0,1x,f5.3,1x,f3.1,']! ')
       read(inc,*,end=30)c0,a0,x0
      endif
      if(command.eq.'q')then
   33  call pgend
       return
      endif
      if(command.eq.'h'.or.command.eq.'p'.or.command.eq.'v')then
       hardcopy=.true.
       goto 34
      endif
      if(command.eq.'?')then
       write(outc,*)'f     Enter new Lorentzian parameters'
       write(outc,*)'h,p,v Hardcopy'
       write(outc,*)'q     Quit'
       goto 31
      endif
      call pgend
      goto 30 
      end

