PROGRAM INVERSEIF

c  General Data Inversion

      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      character endtxt*15,filnam*30,txt*15,name*3,text(2)*55
      character*30 dfilnam(3)
      dimension x(mpar),dat(ndat),aux1(ndat),aux2(ndat),syn(ndat)
      common /xmap/ xxx(mpar),rho(mpar),ind(mpar),inv(mpar),name(mpar)
      common /endpar/ xx(mpar)
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      common /const/ zero,half,one,two
      common /data/ dt,nd,ndf,nd1,nd2,nf1,nf2
      logical nil
      character typ*3,par*76
      data endtxt /'Final Accuracy.'/

      open(4,file='inverseif.out')
      call param(dfilnam,method,nf1,nf2)
      write(6,*) 'Program INVERSEIF - determine deconvolution filter'
      write(4,*) 'Program INVERSEIF - determine deconvolution filter'
      text(1)=' *** CONGRAD ***   Methode of Conjugate Gradients'
      text(2)=' *** MATINV ***   Stabilized Matrix Inversion'
      write(6,*) text(method)
      write(4,*) text(method)
      write(6,*)
      write(4,*) 
      write(6,1) mp,m
      write(4,1) mp,m
    1 format(1x,i3,' Parameters, of which ',i2,' are active')
      write(6,'(a,12i5)') '     k = ',(k,k=1,mp)
      write(4,'(a,12i5)') '     k = ',(k,k=1,mp)
      write(6,'(a,12i5)') ' ind(k)= ',(ind(k),k=1,mp)
      write(4,'(a,12i5)') ' ind(k)= ',(ind(k),k=1,mp)
      write(6,'(a,12i5)') ' inv(k)= ',(inv(k),k=1,m)
      write(4,'(a,12i5)') ' inv(k)= ',(inv(k),k=1,m)
      write(6,3) maxit,qac,finstep
      write(4,3) maxit,qac,finstep
    3 format(' maxit=',i3,'   qac=',f12.9,'   finstep=',f9.6)
      write(4,*)
      write(6,*)

      write(*,*) 'max. number of samples per input file is ',ndat
      call indat(1,dfilnam,dat,aux1,aux2)
      call indat(2,dfilnam,dat,aux1,aux2)
      if(nf1.eq.0) nf1=1
      nf1=max(nf1,1)
      if(nf2.eq.0) nf2=nd
      nf2=min(nf2,nd)
      
      write(6,*)
      write(4,*)
      write(6,*) 'trace '//trim(dfilnam(2))//' to be deconvolved into tr
     &ace '//trim(dfilnam(1))
      write(4,*) 'trace '//trim(dfilnam(2))//' to be deconvolved into tr
     &ace '//trim(dfilnam(1))
      write(6,*) 'from sample ',nf1,' to sample ',nf2
      write(4,*) 'from sample ',nf1,' to sample ',nf2


      qn=zero
      do 6 j=1,nd
    6   qn=qn+dat(j)**2

C  prepare the iterative inversion

      mq=0
      iter=0
      step=one
      axi=step
      noimp=0
      mconj=m
      wq=one
      do 111 k=1,mp
  111 x(k)=zero

c  Start Model

      q=quad(x,dat,aux1,aux2,syn)
      filnam='synt.start'
      txt='Synthetic Start'
      call output(filnam,syn,nd,dt,txt)
      write(6,103) (name(k),k=1,mp)
      write(4,103) (name(k),k=1,mp)
  103 format(/' iter',5x,'rms-err',5(9x,a3)/(17x,5(9x,a3)))
      write(6,104)
      write(4,104)
      write(6,104) iter,dsqrt(q),(xxx(k),k=1,mp)
      write(4,104) iter,dsqrt(q),(xxx(k),k=1,mp)
  104 format(i5,6(1x,f11.6)/(17x,5(1x,f11.6)))
      write(6,105) '+-',(rho(k),k=1,mp)
      write(4,105) '+-',(rho(k),k=1,mp)
  105 format(15x,a2,5(1x,f11.6)/(17x,5(1x,f11.6)))
      write(6,104)
      write(4,104)

C  iterative parameter fitting in diminishing steps between rho and finstep * rho

      write(6,103) (name(inv(k)),k=1,m)
      write(4,103) (name(inv(k)),k=1,m)

      do 2 iter=1,maxit

        qalt=q
        
        if(method.eq.1) then
          call congrd(x,q,gnorm,dat,aux1,aux2,syn)
        else
          call matin(x,q,dat,aux1,aux2,syn)
        endif

        write(6,104) iter,dsqrt(q),(x(k),k=1,m)
        write(4,104) iter,dsqrt(q),(x(k),k=1,m)

c  convergence criterion: within m steps of the iteration, no single
c  step has improved the fit by more that qac or has changed the vector
c  of parameters by more than finstep.

        noimp=noimp+1
        if(dabs(qalt-q).gt.qac.or.axi.gt.finstep) noimp=0
        if(noimp.ge.m.or.axi.lt.finstep/100.d0) goto 5
    2 continue

      endtxt='need more iter.'
      iter=maxit

    5 write(6,104)
      write(4,104)
      write(6,109)
      write(4,109)
  109 format(' Final System Parameters:')
      write(6,103) (name(k),k=1,mp)
      write(4,103) (name(k),k=1,mp)
      write(6,104)
      write(4,104)
      write(6,104) iter,dsqrt(q),(xx(k),k=1,mp)
      write(4,104) iter,dsqrt(q),(xx(k),k=1,mp)
      write(6,108) mq,endtxt
      write(4,108) mq,endtxt
  108 format(/i5,' executions of curve. ',a15)
      write(6,104)
      write(4,104)

      filnam='synt.end'
      txt=   'Synthetic final'
      call output(filnam,syn,nd,dt,txt)

      filnam='rest.end'
      txt=   'targ-syn misfit'
      nf3=nf2-nf1+1
      write(par,*) nf3
      do i=nf1,nf2
        aux2(i-nf1+1)=dat(i)-syn(i)
      enddo          
      call trend(nil,typ,par,aux2,nf3)           
      do 7 j=nf2,nf1,-1
    7 aux2(j)=aux2(j-nf1+1)    
      do 17 j=1,nf1-1
   17 aux2(j)=zero
      do 18 j=nf2+1,nd
   18 aux2(j)=zero

      call output(filnam,aux2,nd,dt,txt)

c produce plot-parameter file. Plot aux data only if they have at least nd samples.
        open(9,file='winplot.par')
        write(9,8) nd
    8   format("0,  5,  24.,  16.,  1,  ",i8,", 0.8")
        write(9,'(a)') dfilnam(1)
        write(9,'(a)') dfilnam(2)
        write(9,'(a)') 'synt.start'
        write(9,'(a)') 'synt.end'
        write(9,'(a)') 'rest.end'
        write(9,'(a)') 'signal simulation with program INVERSEIF'
        close(7)
        write(6,*) 'Plot parameters were saved as file winplot.par '
        
      stop
      end

      double precision function quad(x,dat,aux1,aux2,syn)
      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      dimension x(mpar),dat(ndat),aux1(ndat),aux2(ndat),syn(ndat) 
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      common /const/ zero,half,one,two
      common /data/ dt,nd,ndf,nd1,nd2,nf1,nf2

      logical nil
      character typ*3,par*76

      call curve(x,dat,aux1,aux2,syn)
   
c calculate quadratic error

      nf3=nf2-nf1+1
      write(par,*) nf3
      do i=nf1,nf2
        aux2(i-nf1+1)=dat(i)-syn(i)
      enddo
          
      call trend(nil,typ,par,aux2,nf3)
      
      qua=zero
      do 5 i=1,nf3
    5 qua=qua+aux2(i)**2
      quad=qua/qn
      return
      end

      subroutine matin(x,q,dat,aux1,aux2,syn)
      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      dimension x(mpar),a(ndat,mpar),ata(mpar,mpar),rs(mpar),rsw(mpar) 
      dimension syn0(ndat),dx(mpar),x1(mpar),x2(mpar),atw(mpar,mpar)
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      common /const/ zero,half,one,two
      common /data/ dt,nd,ndf,nd1,nd2,nf1,nf2
      dimension dat(ndat),aux1(ndat),aux2(ndat),syn(ndat)
      diag=qn*qac*nd/mp
      
c  partial derivatives

      do 1 i=1,nd
    1   syn0(i)=syn(i)
      xinc=finstep
      do 2 k=1,mp
        x(k)=x(k)+xinc
        call curve(x,dat,aux1,aux2,syn)
        do 3 i=1,nd
    3     a(i,k)=(syn(i)-syn0(i))/xinc
    2     x(k)=x(k)-xinc

c       do 99 i=1,nd
c  99     write(6,9) (a(i,k),k=1,mp)

c build system of equations

      do 5 k=1,mp
      do 5 j=1,mp
        at=zero
        do 6 i=1,nd
    6     at=at+a(i,k)*a(i,j)
    5 ata(j,k)=at
      do 7 j=1,mp
        rs(j)=zero
        do 7 i=1,nd
    7     rs(j)=rs(j)+a(i,j)*(dat(i)-syn0(i))

c  determine best stabilization weight, wq

      wq1=wq
      q0=q
      do 4 k=1,mp
      do 9 j=1,mp
    9 atw(j,k)=ata(j,k)
      rsw(k)=rs(k)
    4 atw(k,k)=atw(k,k)+diag*wq1
      call gauss(atw,mp,rsw,dx)
      ax1=zero
      do 8 j=1,mp
        ax1=ax1+dx(j)**2
    8   x1(j)=x(j)+dx(j)
      ax1=dsqrt(ax1)
      q1= quad(x1,dat,aux1,aux2,syn)
      if(q1.ge.q) then
        wq2=wq1*8
      else
        wq2=wq1/8
      endif
      do 14 k=1,mp
        do 19 j=1,mp
   19     atw(j,k)=ata(j,k)
        rsw(k)=rs(k)
   14   atw(k,k)=atw(k,k)+diag*wq2
      call gauss(atw,mp,rsw,dx)
      ax2=zero
      do 18 j=1,mp
        ax2=ax2+dx(j)**2
   18   x2(j)=x(j)+dx(j)
      ax2=dsqrt(ax2)
      q2= quad(x2,dat,aux1,aux2,syn)

      if(q1.lt.q0.and.q1.le.q2) then
        wq=wq1
        axi=ax1
        q=q1
        do 21 k=1,mp
   21   x(k)=x1(k)
      else if(q2.lt.q0.and.q2.le.q1) then
        wq=wq2
        axi=ax2
        q=q2
        do 22 k=1,mp
   22   x(k)=x2(k)
      else if(q1.gt.q0+qac.and.q2.gt.q0+qac) then
        axi=one
        wq=wq*8
      else
        axi=zero   
      endif
c     write(6,25) '  q0=',q,'  q1=',q1,'  q2=',q2,'  wq=',wq
   25 format(4(a5,f10.6))
      return
      end


      subroutine gauss(aik,m,rs,f)
c  solve linear equations
      parameter(mpar=12)
      implicit double precision (a-h,o-z)
      dimension aik(mpar,mpar),rs(mpar),f(mpar),h(mpar),imax(mpar)
      do 1401 j=1,m
      aikmax=0.d0
      do 1402 k=1,m
      h(k)=aik(j,k)
      if(dabs(h(k)).le.aikmax) go to 1402
      aikmax=dabs(h(k))
      index=k
 1402 continue
      hh=rs(j)
      do 1403 k=1,m
      q=aik(k,index)/h(index)
      do 1404 l=1,m
 1404 aik(k,l)=aik(k,l)-q*h(l)
 1403 rs(k)=rs(k)-q*hh
      do 1405 k=1,m
 1405 aik(j,k)=h(k)
      rs(j)=hh
 1401 imax(j)=index
      do 1406 j=1,m
      index=imax(j)
 1406 f(index)=rs(j)/aik(j,index)
      return
      end


      subroutine congrd(x,q,gnorm,dat,aux1,aux2,syn)
c  Method of Conjugate Gradients after Fletcher and Reeves (1964)
      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      dimension x(mpar),g(mpar),d(mpar),dd(mpar)
      dimension dat(ndat),aux1(ndat),aux2(ndat),syn(ndat)
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      common /const/ zero,half,one,two

c  Partial Derivatives

      do 1 k=1,m
        x(k)=x(k)+finstep
        qq=quad(x,dat,aux1,aux2,syn)
        x(k)=x(k)-finstep-finstep
        qqq=quad(x,dat,aux1,aux2,syn)
        x(k)=x(k)+finstep
    1   g(k)=(qq-qqq)/two/finstep

c  New direction of descent

      if(mconj.ge.m.or.step.ge.one) then
        write(6,6)
        write(4,6)
    6   format(1x)
        gnorm=0.
        do 2 k=1,m
          gnorm=gnorm+g(k)**2
    2     d(k)=-g(k)
        mconj=1

      else
        ga=gnorm
        if(gnorm.eq.zero) then
          axi=gnorm
          return
        endif
        gnorm=zero
        do 3 k=1,m
    3     gnorm=gnorm+g(k)**2
        beta=gnorm/ga
        do 4 k=1,m
    4     d(k)=-g(k)+beta*d(k)
        mconj=mconj+1
      endif

      dlen=0.
      do 5 k=1,m
    5   dlen=dlen+d(k)**2
      dlen=dsqrt(dlen)
      if(dlen.eq.zero) then
        axi=dlen
        return
      endif

      do 8 k=1,m
    8 dd(k)=d(k)/dlen

c  find minimum in this direction

      call mini(x,q,dd,dat,aux1,aux2,syn)
      return
      end


      subroutine mini(x,q,d,dat,aux1,aux2,syn)

c  x und q entered as start values, returned as improved values
c  step size between 1 and finstep

      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      logical tr
      dimension x(mpar),xl(mpar),xm(mpar),xr(mpar),xx(mpar),d(mpar)
      dimension dat(ndat),aux1(ndat),aux2(ndat),syn(ndat)
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      common /const/ zero,half,one,two
      tr=.false.

      call copy(x,q,xm,qm)

c  check one step left and right
  100 do 1 k=1,m
        xl(k)=xm(k)-step*d(k)
    1   xr(k)=xm(k)+step*d(k)

      ql=quad(xl,dat,aux1,aux2,syn)
      qr=quad(xr,dat,aux1,aux2,syn)

      if(ql.lt.qm.and.qr.lt.qm) then
c  shouldn't happen. In case, reduce step size.
        step=half*step
        write(6,'(a,f6.3)') ' qm maximal, half*step=',step
        write(4,'(a,f6.3)') ' qm maximal, half*step=',step
        if(step.lt.finstep) then
          write(6,'(a)') ' Error: Maximum of q encountered. STOP'
          stop
        else
          goto 100
        endif
      endif

c  find minimum of q

 110  if(ql.lt.qm.and.ql.lt.qr) then
c  reverse so that qr is minimal
        call copy(xl,ql,xx,q0)
        call copy(xr,qr,xl,ql)
        call copy(xx,q0,xr,qr)
        do 9 k=1,m
    9     d(k)=-d(k)
      if(tr) write(6,'(a,f6.3)') ' ql minimal,    -step=',step
      if(tr) write(4,'(a,f6.3)') ' ql minimal,    -step=',step
        goto 110

      else if(qr.lt.ql.and.qr.lt.qm) then
c  check if a minimum has been backeted. If so, use quadratic
c  interpolation. If not, double step size.

        xi=(ql-qr)/(ql-two*qm+qr)/two*step
        axi=dabs(xi)
        if(axi.le.step) then
          step=max(axi,finstep)
          do 6 k=1,m
    6       x(k)=xm(k)+xi*d(k)
          q=quad(x,dat,aux1,aux2,syn)
      if(tr) write(6,'(a,f6.3)') ' qr minimal, qI, step=',step
      if(tr) write(4,'(a,f6.3)') ' qr minimal, qI, step=',step
        else
          step=two*step
          call copy(xr,qr,xm,qm)
          do 111 k=1,m
  111       xr(k)=xm(k)+step*d(k)
          qr=quad(xr,dat,aux1,aux2,syn)
      if(tr) write(6,'(a,f6.3)') ' qr minimal, two*step=',step
      if(tr) write(4,'(a,f6.3)') ' qr minimal, two*step=',step
          goto 110
        endif

      else
      
  112   xi=(ql-qr)/(ql-two*qm+qr)/two*step
        axi=dabs(xi)
        step=min(one,max(xi,finstep))
        do 7 k=1,m
    7     x(k)=xm(k)+xi*d(k)
        q=quad(x,dat,aux1,aux2,syn)
      if(tr) write(6,'(a,f6.3)') ' qm minimal, qI, step=',step
      if(tr) write(4,'(a,f6.3)') ' qm minimal, qI, step=',step
      endif

      return
      end


      subroutine copy(x1,q1,x2,q2)
c  copy x and q from 1 to 2
      parameter (mpar=12)
      implicit double precision (a-h,o-z)
      dimension x1(mpar),x2(mpar)
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      do 1 k=1,m
    1   x2(k)=x1(k)
      q2=q1
      return
      end

      subroutine indat(nf,dfilnam,dat,aux1,aux2)
      parameter(ndat=120000)
      implicit double precision(a-h,o-z)
      character dfilnam(3)*30,form*20,text*50
      common /data/ dt,nd,ndf,nd1,nd2,nf1,nf2
      dimension dat(ndat),aux1(ndat),aux2(ndat)
      write(4,*) 'opening data file ',trim(dfilnam(nf))
      write(6,*) 'opening data file ',trim(dfilnam(nf))
      open(7,file=trim(dfilnam(nf)),status='old')
      read(7,'(a50)') text
      write(4,*) text
      write(6,*) text
  102 read(7,'(a50)') text
      if(text(1:1).eq.'%') goto 102
      
      if(nf.eq.1) then
        read(text,'(i10,a20,f10.6)') nd,form,dt
        nd=min(nd,ndat)
        read(7,form) (dat(i),i=1,nd)
      else if(nf.eq.2) then
        read(text,'(i10,a20)') nd1,form
        nd1=min(nd1,ndat)
        read(7,form) (aux1(i),i=1,nd1)
      endif
         
      close(7)
      write(6,*) i-1,' samples were read'
      return
  100 write(6,*) 'Error at sample # ',i
      stop
  101 write(6,*) 'End of file encountered after sample # ',i-1
      stop
      end

      subroutine falt(x,nx,y,ny,z,nz,off)
      implicit double precision (a-h,o-z)
      dimension x(nx),y(ny),z(nz)
      do 1 i=1,nz
        z(i)=off
        do 1 k=max(1,i+1-nx),min(ny,i)
    1     z(i)=z(i)+x(i-k+1)*y(k)
      return
      end

      subroutine output(filnam,x,n,dt,text)
      implicit double precision (a-h,o-z)
      dimension x(n)
      character filnam*30,text*15,form*15
      form='               '
      xmax=0.
      do 2 j=1,n
    2   xmax=max(xmax,abs(x(j)))
      nvor=log10(max(xmax,2.))+1.
      ndec=max(0,10-nvor)
      write(form,3) ndec
    3   format('(5f13.',i1,')')
      open(8,file=trim(filnam))
      write(8,'(a)') text
      write(8,4) n,form,dt,0.,0.
    4 format(i10,a15,5x,3f10.6)
      write(8,form) (x(i),i=1,n)
      close(8)
      write(6,*) text,' was saved as file ',trim(filnam)
      write(6,*)
      return
      end


      subroutine param(dfilnam,method,nf1,nf2)
      parameter (mpar=12)
      implicit double precision (a-h,o-z)
      character dfilnam(3)*30,name*3
      common /xmap/ xxx(mpar),rho(mpar),ind(mpar),inv(mpar),name(mpar)
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      common /const/ zero,half,one,two
      zero=0.d0
      half=0.5d0
      one=1.d0
      two=2.d0
      open(11,file='inverseif.par',status='old')
      read(11,*)
      do nf=1,2
         read(11,'(a30)') dfilnam(nf)
      enddo
      read(11,*) method
      read(11,*) mp
      read(11,*) maxit
      read(11,*) qac
      qac=qac**2
      read(11,*) finstep
      read(11,*) nf1,nf2
      read(11,*)
      do 1 k=1,mp
    1 read(11,*,err=3,end=3) xxx(k),rho(k),name(k)
      m=0
      do k=1,mp
        if(rho(k).gt.zero) then
          m=m+1
          inv(m)=k
          ind(k)=m
        else
          ind(k)=0
        endif
      enddo
      close(11)
      return
    3 write(*,*) 'Error in file inverseif.par. Check # of parameters.'
      stop    
      end

c  ********************************************     SEIFE subroutines

      subroutine rekfl(x,y,n,f0,f1,f2,g1,g2)
c  perform recursive filtering
      implicit double precision (a-h,o-z)
      dimension x(n),y(n)
      xa=0.d0
      xaa=0.d0
      ya=0.d0
      yaa=0.d0
      do 1 j=1,n
      xn=dble(x(j))
      yn=f0*xn+f1*xa+f2*xaa+g1*ya+g2*yaa
      y(j)=sngl(yn)
      xaa=xa
      xa=xn
      yaa=ya
    1 ya=yn
      return
      end

      subroutine rfk(it,t0,h,t0s,hs,f0,f1,f2,g1,g2)
c  determine coefficients for recursive filter
      implicit double precision (a-h,o-z)
      data zero,one,two,four,eight/0.d0,1.d0,2.d0,4.d0,8.d0/
      if(it.ne.1) goto 10
      f0=one/two/t0
      f1=f0
      f2=zero
      g1=one
      g2=zero
      return
   10 zpi=eight*datan(one)
      eps=zpi/t0
      f2=zero
      g2=zero
      if(it.gt.20) goto 20
      g1=(two-eps)/(two+eps)
      if(it.gt.11) goto 12
      f0=eps/(two+eps)
      f1=f0
      goto 14
   12 if(it.gt.12) goto 13
      f0=two/(two+eps)
      f1=-f0
      goto 14
   13 if(it.gt.13) return
      epss=zpi/t0s
      f0=(epss+two)/(eps+two)
      f1=(epss-two)/(eps+two)
   14 return
   20 epsq=eps*eps
      a=one-eps*h+epsq/four
      b=-two+epsq/two
      c=one+eps*h+epsq/four
      g1=-b/c
      g2=-a/c
      if(it.gt.21) goto 22
      f0=epsq/four/c
      f1=f0+f0
      f2=f0
      goto 25
   22 if(it.gt.22) goto 23
      f0=one/c
      f1=-f0-f0
      f2=f0
      goto 25
   23 if(it.gt.23) goto 24
      epss=zpi/t0s
      epssq=epss*epss
      as=one-epss*hs+epssq/four
      bs=-two+epssq/two
      cs=one+epss*hs+epssq/four
      f0=cs/c
      f1=bs/c
      f2=as/c
      goto 25
   24 if(it.gt.24) return
      f0=eps/two/c
      f1=zero
      f2=-f0
   25 return
      end



      subroutine factor(nil,par,x,n)
c  multiply by a constant factor
      implicit double precision (a-h,o-z)
      dimension x(n)
      character*76 par
      logical nil
      read(par,*,err=99,end=99) f
      do 1 j=1,n
    1 x(j)=f*x(j)
c     write(6,'(" fac  ",f10.3)') f
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for fac'
      stop
      end

      subroutine norm(nil,typ,par,x,n)
c  normalize rms, peak, or peak-to-peak amplitude
      implicit double precision (a-h,o-z)
      dimension x(n)
      character typ*3,par*76
      logical nil
      read(par,*,err=99,end=99) fn
      xrms=0.
      xmin=1e24
      xmax=-xmin
      do 1 j=1,n
        xrms=xrms+x(j)**2
        xmin=min(xmin,x(j))
    1   xmax=max(xmax,x(j))
      xrms=sqrt(xrms/n)
      xp=max(abs(xmin),abs(xmax))
      if(typ.eq.'nox') then
      f=2.*fn/(xmax-xmin)
      off=(xmax+xmin)/2.
        do 2 j=1,n
    2     x(j)=f*(x(j)-off)
      else
        if(typ.eq.'nor') then
          f=fn/xrms
        else
          f=fn/xp
        endif
        do 3 j=1,n
    3     x(j)=f*x(j)
      endif
c     write(6,'(1x,a,a,e12.6)') typ,': normalization factor ',f
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine pad(nil,par,x,n)
c  Nullen ergaenzen
      implicit double precision (a-h,o-z)
      dimension x(n)
      character*76 par
      logical nil
      read(par,*,err=99,end=99) nul
      do 1 j=n+1,nul
    1 x(j)=0.
      n=nul
c     write(6,2) nul
    2 format(" pad: Es wurden Nullen aufgefuellt bis Punkt",i8)
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for pad'
      stop
      end



      subroutine zero(nil,par,x,n)
c  null time series from n1 to n2
      implicit double precision (a-h,o-z)
      dimension x(n)
      character*76 par
      logical nil
      read(par,*,err=99,end=99) n1,n2
      do 1 j=n1,n2
    1 x(j)=0.
c     write(6,'(a,2i10)') ' null ',n1,n2
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for null'
      stop
      end

      subroutine deriv(nil,par,x,n,dt)
c  derivative (non-recursive, non-causal, symmetric-difference)
      implicit double precision (a-h,o-z)
      dimension x(n)
      character*76 par
      logical nil
      read(par,*,err=99,end=99) tau
      if(tau.eq.0.) tau=1.
      twodt=2.*dt/tau
      do 1 j=1,n-2
    1 x(j)=(x(j+2)-x(j))/twodt
      do 2 j=n-1,2,-1
    2 x(j)=x(j-1)
      x(n)=x(n-1)
c     write(6,*) 'dif - derivative by symmetric differences'
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for dif'
      stop
      end

      subroutine mean(nil,par,x,n)
c  remove average
      implicit double precision (a-h,o-z)
      character*76 par
      dimension x(n)
      logical nil
      read(par,*,err=99,end=99) n2
      if(n2.lt.1.or.n2.gt.n) n2=n
      sum=0.
      do 1 j=1,n2
    1 sum=sum+x(j)
      sum=sum/n2
      do 2 j=1,n
    2 x(j)=x(j)-sum
c     write(6,'(" avg  ",e11.3)') sum
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for avg'
      stop
      end


      subroutine trend(nil,typ,par,x,n)
c  remove trend
      implicit double precision (a-h,o-z)
      dimension x(n)
      logical nil
      character typ*3,par*76
        read(par,*,err=99,end=99) n2
        if(n2.le.1.or.n2.gt.n) n2=n
        gn = n2
        alpha = 0.5d0*gn*(gn+1.d0)
        beta = (2.d0*gn+1.d0)*(gn+1.d0)*gn/6.d0
        det = gn*beta-alpha*alpha
        sx = 0.d0
        sjx = 0.d0
        do 1001 j=1,n2
          sx = sx+x(j)
 1001     sjx = sjx+x(j)*j
        a = (sx*beta-sjx*alpha)/det
        b = (sjx*gn-sx*alpha)/det
      do 1002 j=1,n
 1002   x(j) = x(j)-a-b*j
c     write(6,'(" tre  ",2f10.3)') a,b
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for ',typ
      stop
      end
      
      subroutine window(nil,typ,par,x,n,dt,tmin,tsec)
c  time window by sample number
      implicit double precision (a-h,o-z)
      dimension x(n)
      character*76 par
      character*3 typ
      logical nil
      if(typ.eq.'skp'.or.typ.eq.'cos') then
        read(par,*,err=99,end=99) np1
c       write(6,*) typ,np1
      else if(typ.eq.'swi') then
        read(par,*,err=99,end=99) sec1,sec2
        np1=max(nint(sec1/dt)+1,1)
        np2=min(nint(sec2/dt),n)
c       write(6,*) typ,(np1-0.5)*dt,(np2-0.5)*dt
        np1=np1-1
        n=np2-np1
      else
        read(par,*,err=99,end=99) np1,np2
        np1=max(np1,1)
        np2=min(np2,n)
c       write(6,*) typ,np1,np2
        np1=np1-1
        n=np2-np1
      endif
      if(typ.ne.'tap') goto 16
      n=np2
      pih=2.*atan(1.)
      fak=pih/(np2-np1)
      do 14 j=np1,np2
   14 x(j)=x(j)*cos(fak*(j-np1))**2
      goto 15
   16 if(typ.ne.'cos') goto 13
      n1=n+1
      pih=2.*atan(1.)
      fak=pih/np1
      do 17 j=1,np1
      taper=sin(fak*j)**2
      x(j)=x(j)*taper
   17 x(n1-j)=x(n1-j)*taper
      goto 15
   13 if(typ.eq.'win'.or.typ.eq.'swi') then
      do 11 j=1,n
   11 x(j)=x(np1+j)
      endif
      if(typ.eq.'sin'.or.typ.eq.'sis') then
      pi=4.*atan(1.)
      f=pi/(n+1)
      w2=sqrt(2.)
      nex=1
      if(typ.eq.'sis') nex=2
      do 12 j=1,n
   12 x(j)=w2*x(np1+j)*sin(f*j)**nex
      endif
      if(typ.eq.'skp') then
      n=n-np1
      do j=1,n
      x(j)=x(np1+j)
      enddo
      endif
      tsec=60.*tmin+tsec+np1*dt
      tmin=int(tsec/60.)
      tsec=tsec-60.*tmin
   15 nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine polytrend(nil,par,x,n)
c  remove polynomial trend
      parameter(npar=12)
      implicit double precision (a-h,o-z)
      dimension x(n)
      dimension b(npar),c(npar,npar),a(npar)
      logical nil
      character*76 par
      read(par,*,err=99,end=99) m
      m=min(m,npar-1)
      fnh=n/2.
      one=1.d0

      do j=1, m+1
        do k=1, m+1
          c(j,k)=0
          do i=1, n
            c(j,k)=c(j,k)+(dble(i)/fnh-one)**(j+k-2)
          enddo
        enddo
        b(j)=0
        do i=1,n
          b(j)=b(j)+(dble(i)/fnh-one)**(j-1)*x(i)
        enddo
      enddo
      call gauss2(c,m+1,npar,b,a)
c     write(6,100) (j-1,a(j),j=1,m+1)
  100 format(1x,i5,e15.6)
      do i=1,n
        xpol=a(m+1)
        do j=m,1,-1
          xpol=xpol*(dble(i)/fnh-one)+a(j)
        enddo
        x(i)=x(i)-sngl(xpol)
      enddo

c     write(6,'(" pol  ",i5)') m
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for pol'
      stop
      end
      
      subroutine gauss2(aik,m,n,rs,f)
c  solve linear equations
      implicit double precision (a-h,o-z)
      dimension aik(n,n),rs(n),f(n),h(14),imax(13)
      do 1401 j=1,m
      aikmax=0.d0
      do 1402 k=1,m
      h(k)=aik(j,k)
      if(abs(h(k)).le.aikmax) go to 1402
      aikmax=abs(h(k))
      index=k
 1402 continue
      h(m+1)=rs(j)
      do 1403 k=1,m
      q=aik(k,index)/h(index)
      do 1404 l=1,m
 1404 aik(k,l)=aik(k,l)-q*h(l)
 1403 rs(k)=rs(k)-q*h(m+1)
      do 1405 k=1,m
 1405 aik(j,k)=h(k)
      rs(j)=h(m+1)
 1401 imax(j)=index
      do 1406 j=1,m
      index=imax(j)
 1406 f(index)=rs(j)/aik(j,index)
      return
      end     
      
      subroutine reverse(nil,par,x,n)
c  reverse signal in time
      implicit double precision (a-h,o-z)
      logical nil
      dimension x(n)
      character*76 par
      do j=1,n/2
        jj=n+1-j
        xx=x(j)
        x(j)=x(jj)
        x(jj)=xx
      enddo
c     write(6,'(" rev - signal was reversed")')
      nil=.false.
      return
      end
c
      subroutine filter(nil,zeile,x,n,dt,typ,par)
      implicit double precision (a-h,o-z)
      character typ*3,typ1*3,typ2*3,par*76,zeile*80
      dimension x(n)
      logical nil
      typ=zeile(1:3)
      par=zeile(5:80)
c  Subtract value of first sample before high-pass filtration
      if(typ(1:1).eq.'h') then
        x0=x(1)
        do i=1,n
          x(i)=x(i)-x0
        enddo
      endif
c  decode filter type and read parameters
      if(typ.eq.'lp1') then
      read(par,*,err=99,end=99) t0
      it=11
      else if(typ.eq.'hp1') then
      read(par,*,err=99,end=99) t0
      it=12
      else if(typ.eq.'lp2') then
      read(par,*,err=99,end=99) t0,h
      it=21
      else if(typ.eq.'hp2') then
      read(par,*,err=99,end=99) t0,h
      it=22
      else if(typ.eq.'bp2') then
      read(par,*,err=99,end=99) t0,h
      it=24
      else if(typ.eq.'int') then
      t0=1.d0
      it= 1
      else if(typ.eq.'he1'.or.typ.eq.'le1') then
      read(par,*,err=99,end=99) t0s,t0
      it=13
      else if(typ.eq.'he2'.or.typ.eq.'le2') then
        read(par,*,err=99,end=99) t0s,hs,t0,h
      it=23
      else if(typ.eq.'sez') then
        read(par,*,err=99,end=99) t0s,hs
        t0=1.e12
        h=1.
        it=23
      else
      goto 2
      endif
c  calculate filter weights
      call rfk(it,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
      if(typ.eq.'le1') then
        fac=t0s/t0
      f0=f0*fac
      f1=f1*fac
      endif
      if(typ.eq.'le2') then
        fac=(t0s/t0)**2
      f0=f0*fac
      f1=f1*fac
      f2=f2*fac
      endif
c  perform recursive filtration
      call rekfl(x,x,n,f0,f1,f2,g1,g2)
    8 format(1x,a3,2x,5f10.3)
      nil=.false.
c  confirm execution
c     if(it.eq.1) write(6,8) typ,t0
c     if(it.eq.11.or.it.eq.12) write(6,8) typ,t0
c     if(it.eq.13) write(6,9) typ,t0s,t0
    9 format(1x,a3,2x,f10.3,f16.3)
c     if(it.eq.21.or.it.eq.22.or.it.eq.24) write(6,8) typ,t0,h
      if(it.eq.23) then
        if(typ.eq.'sez') then
c         write(6,8) typ,t0s,hs
        else
c         write(6,7) typ,t0s,hs,t0,h
    7 format(1x,a3,2x,2f10.3,f16.3,f10.3)
        endif
      endif
      return

c  Butterworth filters
    2 if(typ.eq.'lpb'.or.typ.eq.'hpb') then
        read(par,*,err=99,end=99) t0,m
        mm=m/2
        if(typ.eq.'lpb') then
          it1=11
          it2=21
          typ1='lp1'
          typ2='lp2'
        else
          it1=12
          it2=22
          typ1='hp1'
          typ2='hp2'
        endif
        if(m.gt.2*mm) then
          call rfk(it1,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
          call rekfl(x,x,n,f0,f1,f2,g1,g2)
c         write(6,8) typ1,t0
        endif
        pih=2.d0*datan(1.d0)
        do 3 j=1,mm
          h=dsin(pih*(2*j-1)/m)
          call rfk(it2,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
          call rekfl(x,x,n,f0,f1,f2,g1,g2)
c         write(6,8) typ2,t0,h
    3   continue
      nil=.false.
      endif
      return
   99 write(6,*) ' insufficient or wrong parameters for ',typ
      stop
      end


      subroutine seife(x,n,dt,zeile)
c   Commands for signal processing:
c
c   avg  n: remove average (determined from first n samples)
c   dif  differentiate the signal
c   fac  f: multiply each sample by f
c   int  integrate the signal
c   lin  remove linear trend through first and last sample
c   nor  a: normalize amplitude to a
c   nop  a: normalize peak amplitude to a
c   nox  a: normalize extrema to +-a
c   nul  n1  n2: nulls the signal between samples n1 and n1 (inclusive)
c   pad  n; pad with zeroes to length n
c   pol  n: remove polynomial trend of degree n
c   rev  reverse signal in time (useful for backward filtering)
c   skp  n: skip n samples
c   tap  n1  n2: signal is cosine- tapered between samples n1 and n2
c        (signal is unchanged between samples 1 and n1)
c   tre  n: remove linear trend (determined from first n samples)
c   twi  tmin1  tsec1  tmin2  tsec2: window defined by time after midnight
c   swi  sec1 sec2: window defined by seconds after first sample
c   win  n1  n2: time window from sample n1 to n2 (n1 becomes first sample)
c
c   recursive (IIR) filters:
c
c   filters are specified by a three-character code for the type
c                            and 1 to 4 parameters.
c     type: one of lp1, hp1, le1, he1, lp2, hp2, bp2, le2, he2, lpb, hpb
c       lp = low pass, hp= high pass, bp = band pass
c             parameters: corner period [,damping if second-order]
c       le = inverse low pass, he = inverse high pass (equalizers)
c             parameters: old corner period, [old damping], new corner
c             period, [new damping]
c         1 or 2: number of poles, i.e. order of denominator polynomial
c         b:      Butterworth, specify period and order
c     Example: lpb 30. 6 is a 30 sec, sixth-order Butterworth low-pass filter.
c Additional type: sez, inverse second-order high-pass filter (restitution)
c   sez t0 h: remove seismo response with free period t0 and damping h
c

      parameter(ndat=120000)
      implicit double precision (a-h,o-z)
      dimension x(ndat)
      character typ*3,zeile*80,par*76
      logical nil
      nil=.true.
      tmin=0.d0
      tsec=0.d0
      typ=zeile(1:3)
      par=zeile(5:80)

      if(typ.eq.'fac') call factor(nil,par,x,n)
      if(typ.eq.'nor'.or.typ.eq.'nop'.or.typ.eq.'nox') 
     &                                    call norm(nil,typ,par,x,n)
      if(typ.eq.'dif') call deriv(nil,par,x,n,dt)
      if(typ.eq.'win'.or.typ.eq.'skp'.or.typ.eq.'tap'.or.typ.eq.'swi') 
     & call window(nil,typ,par,x,n,dt,tmin,tsec)
      if(typ.eq.'nul') call zero(nil,par,x,n)
      if(typ.eq.'avg') call mean(nil,par,x,n)
      if(typ.eq.'tre'.or.typ.eq.'lin') call trend(nil,typ,par,x,n)
      if(typ.eq.'pol') call polytrend(nil,par,x,n)
      if(typ.eq.'rev') call reverse(nil,par,x,n)
      if(typ.eq.'pad') call pad(nil,par,x,n)
      if(nil) call filter(nil,zeile,x,n,dt,typ,par)
c  produce error message if code unrecognized
      if(nil) then
        write(6,3) typ
    3   format(' prozedure ',a3,' is undefined --- stop')
        stop
      endif
      end


      subroutine curve(x,dat,aux1,aux2,syn)
      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      character name*3,zeile*80
      common /aux/ step,finstep,axi,qac,qn,wq,mq,m,mp,maxit,mconj
      common /const/ zero,half,one,two
      common /data/ dt,nd,ndf,nd1,nd2,nf1,nf2
      dimension x(mpar),dat(ndat),aux1(ndat),aux2(ndat),syn(ndat)
      common /xmap/ xxx(mpar),rho(mpar),ind(mpar),inv(mpar),name(mpar)
      common /endpar/ xx(mpar)

c  explanation of variables:
c  nd is the number of samples per file
c  dt is the sampling interval for time series, or frequency interval for spectra
c  dat is the signal (data set) to be fitted
c  aux1 (if read) is an auxiliary series of nd1 samples. 
c  May be used for the input signal that produces dat as the output
c  (this is needed when a system, rather than a signal, is modelled)
c  aux2 (if read) is another auxiliary series of nd2 samples. May be anything.
c  syn will be the synthetic signal that is fitted to dat
c  All signals must have the same sampling interval and the same number of samples.
c  x are normalized parameters, xx actual parameters, xxx initial values.
c  mp is the number of parameters (not more than mpar). m is the number of 
c  active parameters (for which a nonzero search range was specified)
c  Active and passive parameters are mapped with the integer-arrays ind and inv.
c  mq counts how often 'curve' was evaluated.
c  We now determine the actual values of the parameters.

      do 3 i=1,mp
        if (ind(i).gt.0) then
          xx(i)=xxx(i)+rho(i)*x(ind(i))
        else
          xx(i)=xxx(i)
        endif
c        write(6,*) 'xx(',i,') = ',xx(i)
        
    3 continue
    
c  the following code defines the problem. You must replace it with
c  your own code. The auxiliary data sets aux1 and aux2 can be used
c  for any purpose; they are ignored in the rest of the program. 
c  In this version of INVERS, aux1 is the reference (input) signal
c  and dat is the output signal to be fitted.
c  As a result of 'curve', syn must be the synthetic signal
c  calculated with the actual parameters. This signal will then be
c  fitted to the first data set 'dat' by the inversion algorithm.


c      p3=120.
c      p4=0.705
  
      do j=1,nd
        aux2(j)=aux1(j)
      enddo
      
      
      write(zeile,1) 'fac ',xx(1)
c      write(*,*) zeile
      call seife(aux2,nd,dt,zeile)

      write(zeile,1) 'he2 ',xx(2),xx(3),xx(4),xx(5)
c      write(*,*) zeile
      call seife(aux2,nd,dt,zeile)

      write(zeile,1) 'tre ', 0.d0
c      write(*,*) zeile
    1 format(a3,4(1x,e18.12))
      call seife(aux2,nd,dt,zeile)
      
      do j=1,nd
        syn(j)=aux2(j)
      enddo
      

      mq=mq+1
      return
      end