anchor.diff:  !  Splitting operators are applied in order
   
      use FFFTW
  -   use f90sac
   
      implicit none
   
  
positive.code: !     spol:     Source polarisation measured from N.
  !  E,N,Z if not already allocated
  
     use f90sac
  
     implicit none
  
     type(SACtrace), intent(out)   :: E, N, Z
     real, intent(in), optional :: freq, delta, noise, spol
     character(len=*), intent(in), optional :: wavetype
     real(r4), allocatable :: trace(:), random_noise(:)
  
negative.code:    if (present(wavetype)) wavetype_in = wavetype(1:1)
  
     ! Make trace
     call sw_create_wave(E_misfit,N_misfit,Z_misfit, freq=f, delta=delta_in, &
        noise=noise_in, wavetype=wavetype_in, spol=spol)
     call debug_writetrace('wave')
     ! Apply the first splitting operator, then remove the second operator
     call sw_splitN(E_misfit,N_misfit,Z_misfit, 1, (/phi1/), (/ dt1/))
     call debug_writetrace('wave_s1')
     call sw_splitN(E_misfit,N_misfit,Z_misfit, 1, (/phi2/), (/-dt2/))
     call debug_writetrace('wave_s1_-s2')
     ! Find covariance matrix and get misfit as ratio of smaller to larger
     ! eigenvalues of covariance matrix
     call f90sac_covar2(N_misfit, E_misfit, c)
     eig = sw_eig2x2(c)
     misfit = minval(eig)/maxval(eig)
  
     ! Now do the same but the other way round
     call sw_create_wave(E_misfit,N_misfit,Z_misfit, freq=f, delta=delta_in, &
        noise=noise_in, wavetype=wavetype_in, spol=spol)
     call sw_splitN(E_misfit,N_misfit,Z_misfit, 1, (/phi2/), (/ dt2/))
     call debug_writetrace('wave_s2')
     call sw_splitN(E_misfit,N_misfit,Z_misfit, 1, (/phi1/), (/-dt1/))
     call debug_writetrace('wave_s2_-s1')
     call f90sac_covar2(N_misfit, E_misfit, c)
     eig = sw_eig2x2(c)
     misfit = (misfit + minval(eig)/maxval(eig))/2.
  
  contains
  
  