PROGRAM INVERS

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
      data endtxt /'Final Accuracy.'/

      open(4,file='invers.out')
      call param(ndf,dfilnam,method)
      text(1)=' *** CONGRAD ***   Methode of Conjugate Gradients'
      text(2)=' *** MATINV ***   Stabilized Matrix Inversion'
      write(6,*)
      write(6,*) text(method)
      write(4,*) text(method)
 
      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)
      if(ndf.gt.1) call indat(2,dfilnam,dat,aux1,aux2)
      if(ndf.gt.2) call indat(3,dfilnam,dat,aux1,aux2)
      
      filnam='orig.dat'
      txt='original data'
      call output(filnam,dat,nd,dt,txt)
      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=   'data-syn misfit'
      do 7 j=1,nd
    7 dat(j)=dat(j)-syn(j)
      call output(filnam,dat,nd,dt,txt)

c produce plot-parameter file. Plot aux data only if they have at least nd samples.
        open(9,file='winplot.par')
        nplt=ndf+3
        if(ndf.gt.1.and.nd1.lt.nd) nplt=nplt-1
        if(ndf.gt.2.and.nd2.lt.nd) nplt=nplt-1
        write(9,8) nplt,24.,16.,nd
    8   format("0,",i3,",",f5.1,",",f5.1,",","  1,",i8,", 0.8")
        write(9,'(a)') 'orig.dat'
        if(ndf.gt.1.and.nd1.ge.nd) write(9,'(a)') dfilnam(2)
        if(ndf.gt.2.and.nd2.ge.nd) write(9,'(a)') dfilnam(3)
        write(9,'(a)') 'synt.start'
        write(9,'(a)') 'synt.end'
        write(9,'(a)') 'rest.end'
        write(9,'(a)') 'signal simulation with program INVERS'
        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

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

      qua=zero
      do 5 i=1,nd
    5 qua=qua+(syn(i)-dat(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
      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 nach Fletcher und Reeves (1964)
      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      dimension dat(ndat),aux1(ndat),aux2(ndat),syn(ndat)
      dimension x(mpar),g(mpar),d(mpar),dd(mpar)
      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=zero
        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
      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)
      else if(nf.eq.3) then
        read(text,'(i10,a20)') nd2,form
        nd2=min(nd2,ndat)
        read(7,form) (aux2(i),i=1,nd2)
      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='(5f12.3)'
      open(8,file=trim(filnam))
      write(8,'(a)') text
      write(8,3) n,form,dt,0.,0.
    3 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(ndf,dfilnam,method)
      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='invers.par',status='old')
      read(11,*)
      read(11,*) ndf
      do nf=1,ndf
         read(11,'(a30)') dfilnam(nf)
      enddo
      read(11,*,err=2) method
      read(11,*) mp
      read(11,*) maxit
      read(11,*) qac
      qac=qac**2
      read(11,*) finstep
      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
    2 write(*,*) 'Error in file invers.par. Check number of data files.'
      stop
    3 write(*,*) 'Error in file invers.par. Check number of parameters.'
      stop    
      end


      subroutine curve(x,dat,aux1,aux2,syn)
      parameter (mpar=12,ndat=120000)
      implicit double precision (a-h,o-z)
      character*3 name
      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
      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  ndf data files were read
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
    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  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  

c  xx(1)=frequency of line 1
c  xx(2)=amplitude of line 1
c  xx(3)=damping   of line 1
c  xx(4..6) = same for line 2
c  xx(7..9) = same for line 3
c  dt is the frequency increment of the data series

      do j=1,nd
        f=j*dt
        syn(j)=0.d0
        do line=1,3
          f0=xx(3*line-2)
          g=xx(3*line-1)
          h=xx(3*line)
          syn(j)=syn(j)+g**2/(1.d0+(f0/f-f/f0)**2/4.d0/h**2)
        enddo
        
        if(ndf.eq.1) then
          syn(j)=dsqrt(syn(j))
        else
          syn(j)=dsqrt(syn(j))*aux1(j)
        endif
      enddo

      mq=mq+1
      return
      end