      subroutine unclip (x,n,w)
c
c x(n) - real*8 - time series input, ready with gaps
c                 On return, has the gaps filled in
c w(n) - real*8 - The predictions, returned (may have missing values
c                 where the valid data streches are long. 
c
c Limits (work areas fixed size in core / common ):
c 10000 gaps
c  1000 valid samples around each gap.  
c
c Required subroutines
c valrec -  logical function, tests whether the argument is a valid
c           sample of a time series (not the symbolic missing-record value).
c mrs_val - real*8 function, no argument, returns the missing-record symbol. 
c
c Control: 
c subroutine set_dfpmin_tol   (x) - set the tolerance (real*8)
c      entry set_dfpmin_iter  (i) - maximum n.of iterations (integer)
c      entry set_dfpmin_solve (i) - 4 or 5 to solve for the 4 linear
c                                   parameters or, in addition, also
c                                   for the frequency parameter.
c                                   (not recommended, but default (!))
c      entry set_dfpmin_prt   (q) - more diagn. printing (logical)
c      entry set_dfpmin_debug (q) - debug printing (largely unused)
c      entry set_dfpmin_break (q) - Subroutines that are called from
c                                   dfpmin or further down that road
c                                   may call this entry with .true.
c                                   to stop the iterations. 
c
c Defaults:
c Tolerance : 1.0D-6 
c Iter      :     20
c Print     : .false.
c Debug     : .false.
c Nsolve    :      5
c
c 
      implicit real*8 (a-h,o-z)
      parameter (ndim_unclip=1000, nseg=10000)
      logical   valrec,qv,qbreak,qprint,qdbg
      real*8    mrs_val
      dimension x(n),w(n),parc(5)
      dimension xvalf(nseg),xvalb(nseg),mndxf(nseg),mndxb(nseg)
     &,         ibs(nseg),ies(nseg),nval(nseg),miss(nseg)
      common /cunclip/  pi,par(5),m1,m2
      common /csunclip/ nj,indx(ndim_unclip)
     &,                 xs(ndim_unclip),xf(ndim_unclip)
      common /cdfpmin/  gtolc,iterc,nsolvep
      common /cbreak_fpmin/ qbreak,qprint,qdbg

      external chi2_unclip, dchi2_unclip, funclip, dfunclip
      
      pi=4.d0*datan(1.d0)
      qv=valrec(x(1))
      if (qv) then
         is=1
      else
         is=0
      endif
      qbreak=.false.

      do i=2,n
         if (qv.and..not.valrec(x(i))) then
c begin gap
            if (is.gt.0) then
               miss(is)=miss(is)+1
               ies(is)=i-1
            endif
            qv=.false.
         elseif (.not.qv.and.valrec(x(i))) then
c begin segment
            is=is+1
            ibs(is)=i
            if (is.gt.nseg) then
               print *,'<UNCLIP>T> Too many gaps:',is
               stop
            endif
            miss(is)=0
            nval(is)=1
            qv=.true.
         elseif (qv.and.valrec(x(i))) then
            if (is.gt.0) nval(is)=nval(is)+1 
            continue ! have to scan the segments in another way
         elseif (is.gt.0) then
            miss(is)=miss(is)+1
         endif
         w(i)=mrs_val()
      enddo
      w(1)=mrs_val()

      print '(a,i6)',' <UNCLIP>>> Number of segments: ',is

      do ks=1,is-1
         ma=ibs(ks)+1
         mb=ies(ks)
         vmaxchb=0.d0
         vmaxchf=0.d0
         wch=x(mb)-x(mb-1)
         do i=mb,ma,-1
            vch=x(i)-x(i-1)
            if (dabs(vch).gt.vmaxchb) then
               mndxb(ks)=i
               vmaxchb=dabs(vch)
               xvalb(ks)=x(i)
            elseif (dabs(vch).lt.vmaxchb/10.d0
     &          .or.vch*wch.le.0.d0) then
               goto 101
            endif
         enddo
 101     wch=x(ma)-x(ma-1)
         do i=ma,mb
            vch=x(i)-x(i-1)
            if (dabs(vch).gt.vmaxchf) then
               mndxf(ks)=i
               vmaxchf=dabs(vch)
               xvalf(ks)=x(i)
            elseif (dabs(vch).lt.vmaxchf/10.d0
     &          .or.vch*wch.le.0.d0) then
               goto 102
            endif
         enddo
 102     continue
      enddo

      print '(a)',' <UNCLIP>>> (Pivot f|b, miss, val) of the first 100:'
      print '((12h <UNCLIP>>> ,5(2i7,2i4)))'
     &,      (mndxf(i),mndxb(i),miss(i),nval(i),i=1,min0(is,100))

      do ks=1,is-1
c now there are IS segments. Store start at INDX(is), end at JNDX(is) and
c max change at MNDX(is) with value VMAXCH(is) and curve value xval(is)
c
c interpolation formula
c f(i) = A(i) sin (pi * ( i - mndx(ks)) / (mndx(ks+1)-mndx(ks)) 
c      + ( xval(ks+1) * ( i - mndx(ks)) + xval(ks) * ( mndx(ks+1) - i ))
c        / (mndx(ks+1)-mndx(ks))
c where 
c A(i) = ( A2 * ( i - mndx(ks)) + A1 * ( mndx(ks+1) - i )) 
c        / (mndx(ks+1)-mndx(ks))   
c 
c Do nonlin least-square fit to valid records between mndx(ks) and mndx(ks+1)
c with u as a parameter on frequency.
c
c f[i_] := (a2 (i - mndx) + a1 (mndxp - i))/(mndxp - mndx) 
c   Sin [\[Pi] u (i - mndx)/(mndxp - mndx)] + (x2 (i - mndx) + 
c     x1 (mndxp - i))/(mndxp - mndx)
c mndx  -> mndx(ks)
c mndxp -> mndx(ks+1)
c
c Solve for parameters a1,a2,x1,x2,u
c
c df/da1 = ((-i + m2) Sin[((i - m1) \[Pi] u)/(-m1 + m2)])/(-m1 + m2)
c df/da2 =  ((i - m1) Sin[((i - m1) \[Pi] u)/(-m1 + m2)])/(-m1 + m2)
c df/dx1 =  (-i + m2)/(-m1 + m2)
c df/dx2 =   (i - m1)/(-m1 + m2)
c df/du  =  ((i - m1) (a2 (i - m1) + a1 (-i + m2)) \[Pi] 
c          Cos[((i - m1) \[Pi] u)/(-m1 + m2)])/(-m1 + m2)^2
c
         j=0
         do i=mndxb(ks),mndxf(ks+1)
            if (valrec(x(i))) then
               j=j+1
               if (j.gt.ndim_unclip) then
                  print *
     &,'<UNCLIP>T> Not enough space in COMMON /CSUNCLIP/'
                  print *
     &,'<UNCLIP>I> Increase NDIM_UNCLIP in sas/p/unclip.f (3x)'
                  stop
               endif
               xs(j)=x(i)
               xf(j)=funclip(i)
               indx(j)=i
            endif
         enddo
         nj=j
         par(3)=xvalf(ks)
         par(4)=xvalb(ks+1)
         par(1)=0.9d0
         par(2)=1.1d0
         par(5)=1.0d0
         m1=mndxb(ks)
         m2=mndxf(ks+1)
c        
c Starting solution
         do i=1,5
            parc(i)=par(i)
         enddo
         
         np=nsolvep
         gtol=gtolc
         iter=iterc
         qbreak=.false.
         if (qprint) 
     &        print *,'<unclip>>> calling dfpmin, gtol,iter=',gtol,iter
         call dfpmin (parc,np,gtol,iter,fret, chi2_unclip, dchi2_unclip)
         if (qprint) 
     &        print *,'<unclip>>> fret,sol:',fret,par
         do i=1,5
            par(i)=parc(i)
         enddo

         kount=0
         do i=m1-1,m2+1
            w(i)=funclip(i)
            if (.not.valrec(x(i))) then
               x(i)=w(i)
               kount=kount+1
            endif
         enddo
         if (qprint) 
     &        print '(a,i5,i7)'
     &,             ' <unclip>>> segment, replaced:',ks,kount
      enddo
c
c Finally, trailing gap filling by sample-and-hold 
      xl=x(m2)
      do i=m2,n
         if (valrec(x(i))) then
            xl=x(i)
         else
            x(i)=xl
         endif
      enddo
      return
      end

      block data bdfpmin
      real*8 gtolc
      logical qbreak,qprint,qdbg
      common /cdfpmin/ gtolc,iterc,nsolvep
      common /cbreak_fpmin/ qbreak,qprint,qdbg
      data gtol/1.d-6/, iterc/20/, nsolvep/5/
      data qbreak,qprint,qdbg/3*.false./
      end

      subroutine set_dfpmin_tol (x)
      real*8 x,gtolc
      logical q,qbreak,qprint,qdbg
      common /cdfpmin/ gtolc,iterc,nsolvep
      common /cbreak_fpmin/ qbreak,qprint,qdbg

      gtolc=x
      return
      
      entry set_dfpmin_iter (i)
      iterc=i
      return

      entry set_dfpmin_solve (i)
      nsolvep=i
      return

      entry set_dfpmin_prt (q)
      qprint=q
      return

      entry set_dfpmin_debug (q)
      qdbg=q
      return
      
      entry set_dfpmin_break (q)
      qbreak=q
      return
      end

      double precision function funclip (i)
c
c The interpolation formula
c
      implicit real*8 (a-h,o-z)
      dimension par(5)
      common /cunclip/  pi,par,m1,m2
      
      funclip=(par(2)*(i - m1) + par(1)*(m2 - i))/dfloat(m2 - m1) 
     & *dsin(pi * par(5) * (i - m1)/dfloat(m2 - m1)) 
     & + (par(4) * (i - m1) + par(3)* (m2 - i))/dfloat(m2 - m1)
      return
      end

      double precision function dfunclip (k,i)
c
c The derivative of the interpolated value w.r.t. the parameters
c k=1 - a1 (sine amplitude at start of gap)
c   2 - a2 (... at end)
c   3 - x1 (ordinate at first clamping point = ramp start)
c   4 - x2 (... at second ... ramp end)
c   5 - u  (frequency parameter
c
      implicit real*8 (a-h,o-z)
      dimension par(5)
      common /cunclip/  pi,par,m1,m2
      
      if (k.eq.1) then
         dfunclip=((-i + m2)*dsin(((i - m1)*pi*par(5))
     &            /dfloat(-m1 + m2)))/dfloat(-m1 + m2)
      elseif (k.eq.2) then
         dfunclip=(( i - m1)*dsin(((i - m1)*pi*par(5))
     &            /dfloat(-m1 + m2)))/dfloat(-m1 + m2)
      elseif (k.eq.3) then
         dfunclip=  dfloat(-i + m2)/dfloat(-m1 + m2)
      elseif (k.eq.4) then
         dfunclip=  dfloat( i - m1)/dfloat(-m1 + m2)
      elseif (k.eq.5) then
         dfunclip=(dfloat(i - m1)*( par(2)*(i - m1) + par(1)*(-i + m2))
     &            *pi*dcos(((i - m1)*pi*par(5))/dfloat(-m1 + m2)))
     &             /dfloat(-m1 + m2)**2
      endif
      return
      end

      double precision function chi2_unclip (parc,np)
      implicit real*8 (a-h,o-z)
      parameter (ndim_unclip=1000, nseg=10000)
      dimension parc(5)
      common /cunclip/  pi,par(5),m1,m2
      common /csunclip/ nj,indx(ndim_unclip)
     &,                 xs(ndim_unclip),xf(ndim_unclip)

      chi2=0.d0
      do i=1,np
         par(i)=parc(i)
      enddo

      do j=1,nj
         i=indx(j)
         xf(j)=funclip(i)
         chi2=chi2+(xf(j)-xs(j))**2
      enddo
      chi2_unclip = chi2/nj
      return
      end

      subroutine dchi2_unclip (parc,np,df)
      implicit real*8 (a-h,o-z)
      parameter (ndim_unclip=1000, nseg=10000)
      dimension parc(np),df(np)
      common /cunclip/  pi,par(5),m1,m2
      common /csunclip/ nj,indx(ndim_unclip)
     &,                 xs(ndim_unclip),xf(ndim_unclip)

      do i=1,np
         df(i)=0.d0
         par(i)=parc(i)
      enddo

      do j=1,nj
         i=indx(j)
         xf(j)=funclip(i)
         do k=1,np
            df(k)=df(k)+2.d0*(xf(j)-xs(j))*dfunclip(k,i)
         enddo
      enddo
      
      do i=1,np
         df(i)=df(i)/nj
      enddo
      return
      end

         
      
