      subroutine huberwt(nsect,scale,lwt,res,w)
c
c  birrp
c
c  computes huber weights based on regression residuals
c
c  input variables
c    nsect-number of sections
c    scale-robust estimate of scale
c    lwt-logical vector, .true. if datum is used
c    res-vector of absolute value of residuals
c  output variable
c    w-vector of weights
c
      logical lwt
      dimension lwt(*),res(*),w(*)
      do 10 i=1,nsect
        if(.not.lwt(i))then
          w(i)=0.
        else
          t=res(i)/scale
          if(t.le.1.5)then
            w(i)=1.
          else
            w(i)=1.5/t
          endif
        endif
   10   continue
      return
      end
      subroutine expwt(nsect,u0,scale,lwt,res,w)
c
c  birrp
c
c  computes exponential (thomson) weights based on regression residuals
c
c  input variables
c    nsct-number of sections
c    u0-parameter for exponential weights
c    scale-robust estimate of scale
c    lwt-logical vector, .true. if datum is used
c    res-vector of absolute value of residuals
c  output variable
c    w-vector of weights
      logical lwt
      dimension lwt(*),res(*),w(*)
      small=r1mach(1)
      big=log(-log(small)/.05129329)
c  coefficient sets the 0.95 point to res/scale=u0
      sw=exp(-.05129329*exp(-u0**5))
      do 10 i=1,nsect
        if(.not.lwt(i))then
          w(i)=0.
        else
          t=u0**4*(res(i)/scale-u0)
          if(t.lt.big)then
            w(i)=exp(-.05129329*exp(t))/sw
          else
            w(i)=0.
          endif
        endif
   10   continue
      return
      end
      subroutine hatwt(nsect,ninp,ainl,ainu,llo,lhi,hat,what)
c
c  birrp
c
c  computes leverage weights based on size of hat matrix diagonal relative
c  to its expected value
c
c  input variables
c    nsect-number of sections
c    ninp-number of inputs
c    ainl-lower cutoff
c    ainu-upper cutoff
c    llo-logical vector denoting data deleted at lower end of distribution,
c        changed on return
c    lhi-logical vector denoting data deleted at upper end of distribution,
c        changed on return
c    hat-vector of hat matrix diagonal entries
c  output variables
c    what-vector of weights
c
      logical llo,lhi
      dimension llo(*),lhi(*)
      dimension hat(*),what(*)
      small=r1mach(1)
      big=log(-log(small)/.05129329)
      nw=0
      do 5 i=1,nsect
        if(what(i).ne.0.)nw=nw+1
    5   continue
      xnw=nw
      if(ainu.le.0.)goto 20
      ains=ainu
      if(ains.lt.1.)ains=1./ains
      swhat=exp(-.05129329*exp(-ainu*ains**4))
      do 10 i=1,nsect
        if(what(i).eq.0.)goto 10
        t=ains**4*(xnw*hat(i)/ninp-ainu)
        if(t.lt.big)then
          owhat=what(i)
          what(i)=what(i)*exp(-.05129329*exp(t))/swhat
          if((owhat.ne.0.).and.(what(i).eq.0.))lhi(i)=.true.
        else
          what(i)=0.
          lhi(i)=.true.
        endif
  10    continue
  20  if(ainl.le.0.)return
      ains=ainl
      if(ains.gt.1.)ains=1./ains
      swhat=exp(-.05129329*exp(-ainl/ains**4))
      do 30 i=1,nsect
        if(what(i).eq.0.)goto 30
        t=(xnw*hat(i)/ninp-ainl)/ains**4
        if(t.lt.big)then
          owhat=what(i)
          what(i)=what(i)*(1.-exp(-.05129329*exp(t))/swhat)
          if((owhat.ne.0.).and.(what(i).eq.0.))llo(i)=.true.
        endif
  30    continue
      return
      end
      subroutine mad(nsect,res,scale)
c
c  birrp
c
c  computes median absolute deviation for a Rayleigh variate
c
c  input variables:
c
c  nsect--number of entries in res
c  res--vector of residual magnitudes (overwritten on return)
c
c  output variables:
c
c  scale-mad of entries in res
c
      dimension res(*)
      call ssort(res,res,nsect,1)
      ni = nsect
      do 20 i=1,nsect
        if(res(i).gt.0.)then
          ni=i
          goto 21
        endif
   20   continue
   21 nl=nsect-ni+1
      if(mod(nl,2).eq.0.)then
        xmed=(res(nl/2+ni-1)+res(nl/2+ni))/2.
      else
        xmed=res(nl/2+ni)
      endif
      do 30 i=ni,nsect
   30   res(i)=abs(res(i)-xmed)
      call ssort(res(ni),res(ni),nl,1)
      if(mod(nl,2).eq.0)then
        scale=(res(nl/2+ni-1)+res(nl/2+ni))/.8969
      else
        scale=res(nl/2+ni)/.44845
      endif
      return
      end
      subroutine hatdiag(nsct,ninp,ftsec,wt,hat,info)
c
c  birrp
c
c  computes diagonal of hat matrix 
c
c  input variables:
c
c  nsct--number of sections
c  ninp--number of variables
c  ftsec--complex array containing nsct section ffts for ninp variables
c  wt--weight vector
c
c  output variables:
c
c  hat--vector containing hat matrix diagonal
c  info--error flag, 0 for normal, otherwise as returned by zhifa
c
      include 'parameters.h'
      double precision det
      complex ftsec
      double complex a,atemp,zz
      dimension ftsec(nsectm,*),wt(*),hat(*)
      dimension a(nrefm,nrefm),ipvt(nrefm),det(2),zz(nrefm)
      info=0
      do 10 i=1,ninp
        do 10 j=i,ninp
          a(i,j)=0.
          do 10 ii=1,nsct
   10       a(i,j)=a(i,j)+wt(ii)**2*conjg(ftsec(ii,i))*ftsec(ii,j)
      call zhifa(a,nrefm,ninp,ipvt,info)
      job=001
      call zhidi(a,nrefm,ninp,ipvt,det,inert,zz,job)
      do 20 i=1,ninp
        do 20 j=i+1,ninp
   20     a(j,i)=conjg(a(i,j))
      do 30 j=1,nsct
        hat(j)=0.
        do 30 i=1,ninp
          atemp=0.
          do 25 jj=1,ninp
   25       atemp=atemp+a(i,jj)*wt(j)*conjg(ftsec(j,jj))
   30     hat(j)=hat(j)+wt(j)*real(ftsec(j,i)*atemp)
      return
      end
      subroutine residual(nsect,nout,ninp,nn,ftsec,h,wt,nw,res,resq)
c
c  birrp
c
c  returns residuals and residual sum of squares
c
c  input variables
c
c   nsect--number of sections
c   nout--number of output variables
c   ninp--number of input variables
c   nn--index of output variable
c   ftsec--complex nsect by nout+ninp array containing fourier transforms
c   h--double complex vector containing ninp transfer function estimates
c   wt--vector of weights
c
c  returned variables
c
c  nw--number of residuals with nonzero weights
c  res--vector of absolute residuals
c  resq--normalized residual sum of squares
c
      include 'parameters.h'
      complex ftsec,resq
      double complex h,ctmp,resq1
      dimension ftsec(nsectm,*),h(*),wt(*),res(*)
      resq1=0.
      nw=0
      do 20 n=1,nsect
        ctmp=ftsec(n,nn)
        if(wt(n).gt.0.)nw=nw+1
        do 10 i=1,ninp
   10     ctmp=ctmp-h(i)*ftsec(n,nout+i)
        resq1=resq1+wt(n)**2*dcmplx(dreal(ctmp)**2,dimag(ctmp)**2)
   20   res(n)=abs(ctmp)
      resq=resq1/nw
      return
      end
      subroutine delone1(nsect,ninp,nref,nn,ftsec,wt,aa,bb,hd1,info)
c
c  birrp
c
c  computes delete one magnetic transfer transfer functions
c
c  input variables:
c
c  nsect--number of sections
c  ninp-number of local magnetic field time series
c  nref--number of reference time series
c  nn--index of local time series
c  ftsec--complex array of section ffts
c  wt--vector of weights
c  aa-double complex array containing input spectral matrix
c  bb-double complex vector containing input/output spectral matrix
c
c  output variables:
c
c  hd1--complex array of delete one responses
c  info-error flag, 0 for normal return, as returned by zgeco otherwise
c
      include 'parameters.h'
      complex ftsec,hd1
      double complex aa,bb,a,b,zz
      double precision rcond
      dimension ftsec(nsectm,*),aa(nrefm,*),bb(*),
     $          wt(*),hd1(nrefm+1,*)
      dimension a(nrefm,nrefm),b(nrefm),ipvt(nrefm),zz(nrefm)
      info=0
      do 10 ii=1,nsect
        do 5 i=1,nref
          b(i)=bb(i)-wt(ii)**2*conjg(ftsec(ii,ninp+i))*
     $         ftsec(ii,nn)
          do 5 j=1,nref
    5       a(i,j)=aa(i,j)-wt(ii)**2*conjg(ftsec(ii,ninp+i))*
     $             ftsec(ii,ninp+j)
          call zgeco(a,nrefm,nref,ipvt,rcond,zz)
          if(1.d0+rcond.eq.1.d0)info=1
          job=0
          call zgesl(a,nrefm,nref,ipvt,b,job)
          do 10 i=1,nref
  10        hd1(i,ii)=b(i)
      return
      end
      subroutine delone2(nsect,nout,ninp,nn,ftsec,wt,aa,bb,hd1,info)
c
c  birrp
c
c  computes delete one transfer functions
c
c  input variables:
c
c  nsect--number of sections
c  nout-number of output time series
c  ninp--number of input time series
c  nn--index of output time series
c  ftsec--complex array of nsect section ffts
c  wt--vector of weights
c  aa-double complex array containing input spectral matrix
c  bb-double complex vector containing input/output spectral matrix
c
c  output variables:
c
c  hd1--complex array of delete one responses
c  info-error flag, 0 for normal return, as returned by zgeco otherwise
c
      include 'parameters.h'
      complex ftsec,hd1
      double complex aa,bb,a,b,zz
      double precision rcond
      dimension ftsec(nsectm,*),aa(nrefm,*),bb(*),
     $          wt(*),hd1(nrefm+1,*)
      dimension a(ninpm,ninpm),b(ninpm),ipvt(ninpm),zz(ninpm)
      info=0
      do 10 ii=1,nsect
        do 5 i=1,ninp
          b(i)=bb(i)-wt(ii)**2*conjg(ftsec(ii,nout+ninp+i))*
     $         ftsec(ii,nn)
          do 5 j=1,ninp
    5       a(i,j)=aa(i,j)-wt(ii)**2*conjg(ftsec(ii,nout+ninp+i))*
     $             ftsec(ii,nout+j)
           call zgeco(a,ninpm,ninp,ipvt,rcond,zz)
          if(1.d0+rcond.eq.1.d0)info=1
          job=0
          call zgesl(a,ninpm,ninp,ipvt,b,job)
          do 10 i=1,ninp
  10        hd1(i,ii)=b(i)
      return
      end