      subroutine qqout(filnam,period,nsect,nsct,nw,ninp,llo,lhi,
     $                 wt,res,hat)
      include 'parameters.h'
      character*(*) filnam
      logical llo,lhi
      dimension llo(*),lhi(*),wt(*),res(*),hat(*)
      dimension temp(nsectm),xindex(nsectm),qq(nsectm),yindex(nsectm)
      open (unit=40,file=filnam,status='unknown')
      xnsect=nsect
      xnw=nw
      m1=0
      m2=0
      do 10 i=1,nsect
        if(llo(i))m1=m1+1
   10   if(lhi(i))m2=m2+1
c  apportion data thresholded or eliminated as outliers between upper and lower ends of
c  beta distribution
      nray=nsect-nw-m1-m2
      nloh=nray/2
      nhih=nray-nloh
c apportion data thresholded or eliminated as leverage points between upper and lower
c  ends of rayleigh distribution
      nbeta=nsect-nsct+m1+m2
      nlor=nbeta/2
      nhir=nbeta-nlor
      nr=nsct-nw-m1-m2
      var=0.
      do 20 i=1,nsect
        temp(i)=(wt(i)*res(i))**2
   20   var=var+temp(i)
c  scale regression residuals
c  set second moment to 2 as expected for rayleigh-distributed variate
      var=sqrt(2.*xnsect/var)
      do 30 i=1,nsect
        xindex(i)=i
        qq(i)=var*sqrt(temp(i))
	yindex(i)=i
        temp(i)=xnsect*hat(i)/ninp
   30   continue
      call ssort(qq,xindex,nsect,2)
      call ssort(temp,yindex,nsect,2)
      do 40 i=1,nsect
        if((qq(i).gt.0.).and.(temp(i).gt.0.))then
          nq=i
          goto 41
        endif
   40   continue
   41 continue
      ft1=0.
      if(nlor.gt.0)then
        t1=sqrt(-2.*log(1.-(nlor-.5)/xnsect))
        ft1=1.-exp(-t1*t1/2.)
      endif
      ft2=1.
      if((nr.gt.0).or.(nhir.gt.0))then
        t2=sqrt(-2.*log(1.-(xnsect-nhir-nr-.5)/xnsect))
        ft2=1.-exp(-t2*t2/2.)-ft1
      endif
      ft3=0.
      if((m1.gt.0).or.(nloh.gt.0))then
        t3=qbeta((float(m1+nloh)-.5)/xnsect,ninp,nsect-ninp)
        ft3=betai(t3,float(ninp),xnsect-float(ninp))
      endif
      ft4=1.
      if((m2.gt.0).or.(nhih.gt.0))then
        t4=qbeta((float(nsect-m2-nhih)-.5)/xnsect,ninp,nsect-ninp)
        ft4=betai(t4,float(ninp),xnsect-float(ninp))-ft3
      endif
      write(40,999)period,0.,0.,0.,0.,0.
      do 50 i=nq,ifix(0.1*xnw)+nq
        p=(float(i-nq+1)-.5)/xnw
   50   write(40,999)sqrt(-2.*log(1.-ft2*p-ft1)),
     $               qq(i),xindex(i),
     $               xnsect*qbeta(ft4*p+ft3,ninp,nsect-ninp)/ninp,
     $               temp(i),yindex(i)
      do 60 i=ifix(0.1*xnw)+nq+1,ifix(0.9*xnw)+nq,
     $        max(nint(log10(xnw)**1.5),1)
        p=(float(i-nq+1)-.5)/xnw
   60   write(40,999)sqrt(-2.*log(1.-ft2*p-ft1)),
     $               qq(i),xindex(i),
     $               xnsect*qbeta(ft4*p+ft3,ninp,nsect-ninp)/ninp,
     $               temp(i),yindex(i)
      do 70 i=ifix(0.9*xnw)+nq+1,nsect
        p=(float(i-nq+1)-.5)/xnw
   70   write(40,999)sqrt(-2.*log(1.-ft2*p-ft1)),
     $               qq(i),xindex(i),
     $               xnsect*qbeta(ft4*p+ft3,ninp,nsect-ninp)/ninp,
     $               temp(i),yindex(i)
      close (unit=40)
  999 format(1x,6g15.7)
      return
      end
      subroutine diagout(filnam,period,nsect,ninp,
     $                   res,pwro,pwri,wt,hat)
      character*(*) filnam
      dimension res(*),pwro(*),pwri(*),wt(*),
     $          hat(*)
      xnsect=nsect
      open (unit=40,file=filnam,status='unknown')
      respm=0.
      pwrom=0.
      pwrim=0.
      do 10 i=1,nsect
        respm=max(respm,abs(res(i))**2)
        pwrom=max(pwrom,pwro(i))
   10   pwrim=max(pwrim,pwri(i))
      write(40,999)period,respm,pwrom,pwrim
      do 20 i=1,nsect
        if(wt(i).gt.0.)write(40,999)xnsect*hat(i)/ninp,
     $                 abs(res(i))**2/respm,pwro(i)/pwrom,
     $                 pwri(i)/pwrim
   20   continue
      close(unit=40)
  999 format(1x,4g15.7)
      return
      end
