      subroutine mse(n,a,x,z,err,mdim,edim,iwin)
c      subroutine msx(ipltall,ierr,ipen,nn1,nn2,iwin,n,a,x,z,eps,
c     1               einc,edec)
c
c  Multi-source Euler deconvolution
c
c  Programmed by Jeff Phillips - October, 1994. (Incomplete)
c
      common /work1/xc(512)
      common /work2/xd(512)
      common /work3/aa(84),uu(84),dum(688)
      common /work4/vv(441),s(21),work(63)
      common /model/nbod,ncor,iopn,delx,ibod
      common /pscale/ xmin,xmax,zmin,zmax
      common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      dimension a(n),x(n),z(n),p(63)
      complex xc,xd,w,aa,uu,vv
      complex b(21),ww(21),wc(21),ap(22),num,den,alpha
      complex wsav(512)
      dimension aksav(512),drsav(512),zavesav(512)
      character sym*1
        data d2r/1.745329e-2/
      index(i,j,ny) = ((i - 1) * ny) + j
c
c      read(5,*) mm,cc
      ierr=0
      nsav=0
      call setfft(1,lc,n,a)
      k=lc/2
      do 5 j=2,k
      xc(j)=xc(j)+xc(j)
    5 xc(lc+2-j)=0.
      call fork(lc,xc,-1.)
c
c  xc contains the analytic signal
c
      call csplinec(lc,delx)
c
c  xd contains the first derivative of the analytic signal
c
      size=.08
      if(err.eq.0.) then
        std=1.e38
      else
        std=100./err
      endif
      if(mdim.eq.3) then
        write(20,210) iwin,std
      else
        write(20,211) iwin,std
      endif
  210 format(i10,' point window, ',f5.1,
     1'% maximum error, floating dimension')
  211 format(i10,' point window, ',f5.1,
     1'% maximum error, fixed dimension')
      write(20,206)
  206 format('         x              z             dim            %err'
     1)
      do 900 ii=1,n-iwin
      xmin1=x(ii)
      xmax1=x(ii+iwin)
      if(xmax1.lt.xmin.or.xmin1.gt.xmax) go to 900
c      do 800 nn=nn1,nn2
c      if(nn.gt.10) go to 990
c      if(nn.lt.1) go to 990
c      noo=0
      zmax1=-1.e38
      zmin1=1.e38
      xcmax=-1.e38
      do 6 i=1,n
      if(x(i).lt.xmin1) go to 6
      if(x(i).gt.xmax1) go to 7
c      noo=noo+1
      if(z(i).gt.zmax1) zmax1=z(i)
      if(z(i).lt.zmin1) zmin1=z(i)
      if(cabs(xc(i)).gt.xcmax) xcmax=cabs(xc(i))
   6  continue
   7  np=mdim+1
      if((np)*iwin.gt.84) go to 990
      if(iwin*iwin.gt.441) go to 992
      xave=(xmax1+xmin1)/2.
      zave=(zmax1+zmin1)/2.
      xran=(xmax1-xmin1)/2.
      zran=(zmax1-zmin1)/2.
      if(zran.gt.xran) xran=zran
      no=0
      do 20 i=1,n
      if(x(i).lt.xmin1) go to 20
      if(x(i).gt.xmax1) go to 30
      no=no+1
      w=cmplx((x(i)-xave),(z(i)-zave))/xran
      do 10 k=1,iwin
      ll=index(k,no,iwin)
      aa(ll)=xd(i)/xcmax
      if(mdim.eq.2) then
        aa(ll+iwin)=1.
        aa(ll+2*iwin)=edim*xc(i)+w*xd(i)
      else
        aa(ll+iwin)=-xc(i)/xcmax
        aa(ll+2*iwin)=1.
        aa(ll+3*iwin)=w*xd(i)
      endif
   10 continue
   20 continue
   30 if(no.lt.(np-1)) go to 991
      call csvd(aa,no,np,s,uu,vv,work)
      call vsb(aa,no,np,s,vv)

c800   continue
900   continue

c******************************************************************************
      subroutine cspline(n,dx)
c
c  returns the first derivative, aa, of the analytic signal xc
c
      common /work1/xc(512)
      common /work2/aa(512)
      common /work3/c(512),dum(512)
      common /work4/d(512)
      complex xc,aa,d
      c(1)=.5
      d(1)=1.5*(xc(2)-xc(1))/dx
      n1=n-1
      do 10 i=2,n1
      dnm=(4.0-c(i-1))*dx
c      print *,'dnm = ',dnm,' d(i-1) = ',d(i-1)
      c(i)=dx/dnm
   10 d(i)=(3.0*(xc(i+1)-xc(i-1)) - dx*d(i-1))/dnm
c      print *,'c(n1) = ',c(n1)
      aa(n)=(3.*(xc(n)-xc(n1))/dx - d(n1))/(2.-c(n1))
      do 9 i=1,n1
      k=n-i
    9 aa(k)=d(k)-c(k)*aa(k+1)
      return
      end
c******************************************************************************
        subroutine csvd(a,m,n,s,u,v,wrk)
        dimension s(n-1),wrk(1)
        complex a(m,n),u(m,m),v(n-1,n-1)
        n2=n
        n3=2*n-1
        call csvd1(a,m,n-1,m,n-1,1,m,n-1,s,u,v,wrk,wrk(n2),wrk(n3))
	return
	end
c******************************************************************************
        subroutine csvd1(a,mmax,nmax,m,n,p,nu,nv,s,u,v,b,c,t)
c  singular value decomposition of a matrix
c  Peter A. Businger, Bell Telephone Laboratories
c  Gene H. Golub, Stanford University
c  Algorithm 358, Collected algorithms from ACM
        complex a(mmax,1), u(mmax,1), v(nmax,1)
        dimension s(n), b(n), c(n), t(n)
	integer p
        complex q,r
	data eta,tol/1.5e-8, 1.e-31/
	np=n+p
        n1=n+1
c        ierr=0
c
c  householder reduction
	c(1)=0.0
	k=1
10	k1=k+1
c
c  elimination of a(i,k), i=k+1...m
	z=0.0
	do 20 i=k,m
20      z = z + real(a(i,k))**2+aimag(a(i,k))**2
	b(k)=0.0
	if(z.le.tol) go to 70
	z=sqrt(z)
	b(k)=z
        w=cabs(a(k,k))
        q=(1.0,0.0)
	if(w.ne.0.0) q=a(k,k)/w
	a(k,k)=q*(z+w)
	if(k.eq.np) go to 70
	do 50 j=k1,np
        q=(0.0,0.0)
	do 30 i=k,m
30      q=q+conjg(a(i,k))*a(i,j)
	q=q/(z*(z+w))
	do 40 i=k,m
40      a(i,j)=a(i,j)-q*a(i,k)
50      continue
c
c  phase transformation
        q=-conjg(a(k,k))/cabs(a(k,k))
	do 60 j=k1,np
60      a(k,j)=q*a(k,j)
c
c  elimination of a(k,j), j=k+2...n
70	if(k.eq.n) go to 140
	z=0.0
	do 80 j=k1,n
80      z = z + real(a(k,j))**2+aimag(a(k,j))**2
	c(k1)=0.0
	if(z.le.tol) go to 130
	z=sqrt(z)
	c(k1)=z
        w=cabs(a(k,k1))
        q=(1.0,0.0)
	if(w.ne.0.0) q=a(k,k1)/w
	a(k,k1)=q*(z+w)
	do 110 i=k1,m
        q=(0.0,0.0)
	do 90 j=k1,n
90      q = q + conjg(a(k,j))*a(i,j)
	q=q/(z*(z+w))
	do 100 j=k1,n
100     a(i,j)=a(i,j)-q*a(k,j)
110	continue
c
c  phase transformation
        q=-conjg(a(k,k1))/cabs(a(k,k1))
	do 120 i=k1,m
120     a(i,k1)=a(i,k1)*q
130	k=k1
	go to 10
c
c  tolerance for negligible elements
140	eps=0.0
	do 150 k=1,n
	s(k)=b(k)
	t(k)=c(k)
150	eps=amax1(eps,s(k)+t(k))
	eps=eps*eta
c
c  initialization of u and v
        if(nu.eq.0) go to 180
        do 170 j=1,nu
	do 160 i=1,m
160     u(i,j)=(0.0,0.0)
170     u(j,j)=(1.0,0.0)
180     if(nv.eq.0) go to 210
	do 200 j=1,nv
	do 190 i=1,n
190     v(i,j)=(0.0,0.0)
200     v(j,j)=(1.0,0.0)
c
c  qr diagonalization
210	do 380 kk=1,n
	k=n1-kk
c
c  test for split
220	do 230 mm=1,k
	m2=k+1-mm
	if(abs(t(m2)).le.eps) go to 290
	if(abs(s(m2-1)).le.eps) go to 240
230	continue
c
c  cancellation of e(m2)
240	cs=0.0
	sn=1.0
	m1=m2-1
	do 280 i=m2,k
	f=sn*t(i)
	t(i)=cs*t(i)
	if(abs(f).le.eps) go to 290
	h=s(i)
	w=sqrt(f*f+h*h)
	s(i)=w
	cs=h/w
	sn=-f/w
	if(nu.eq.0) go to 260
	do 250 j=1,n
        x=real(u(j,m1))
        y=real(u(j,i))
        u(j,m1) = cmplx(x*cs + y*sn,0.0)
250     u(j,i)  = cmplx(y*cs - x*sn,0.0)
260	if(np.eq.n) go to 280
	do 270 j=n1,np
	q=a(m1,j)
	r=a(i,j)
	a(m1,j) = q*cs + r*sn
        a(i,j)  = r*cs - q*sn
270     continue
280	continue
c
c  test for convergence
290	w=s(k)
	if(m2.eq.k) go to 360
c
c  origin shift
	x=s(m2)
	y=s(k-1)
	g=t(k-1)
	h=t(k)
	f=((y-w)*(y+w)+(g-h)*(g+h))/(2.0*h*y)
	g=sqrt(f*f+1.0)
	if(f.lt.0.0) g=-g
	f=((x-w)*(x+w)+(y/(f+g)-h)*h)/x
c
c  qr step
	cs=1.0
	sn=1.0
	m1=m2+1
	do 350 i=m1,k
	g=t(i)
	y=s(i)
	h=sn*g
	g=cs*g
	w=sqrt(h*h+f*f)
	t(i-1)=w
	cs=f/w
	sn=h/w
	f=x*cs+g*sn
	g=g*cs-x*sn
	h=y*sn
	y=y*cs
	if(nv.eq.0) go to 310
	do 300 j=1,n
        x=real(v(j,i-1))
        w=real(v(j,i))
        v(j,i-1)= cmplx(x*cs + w*sn,0.0)
300     v(j,i)  = cmplx(w*cs - x*sn,0.0)
310	w=sqrt(h*h+f*f)
	s(i-1)=w
	cs=f/w
	sn=h/w
	f=cs*g+sn*y
	x=cs*y-sn*g
	if(nu.eq.0) go to 330
	do 320 j=1,n
        y=real(u(j,i-1))
        w=real(u(j,i))
        u(j,i-1) = cmplx(y*cs + w*sn,0.0)
320     u(j,i)   = cmplx(w*cs - y*sn,0.0)
330	if(n.eq.np) go to 350
	do 340 j=n1,np
	q=a(i-1,j)
	r=a(i,j)
	a(i-1,j)=q*cs+r*sn
        a(i,j) = r*cs-q*sn
340     continue
350	continue
	t(m2)=0.0
	t(k)=f
	s(k)=x
	go to 220
c
c  convergence
360	if(w.ge.0.0) go to 380
	s(k)=-w
	if(nv.eq.0) go to 380
	do 370 j=1,n
370	v(j,k)=-v(j,k)
380	continue
c
c  sort singular values
	do 450 k=1,n
	g=-1.0
	j=k
	do 390 i=k,n
	if(s(i).le.g) go to 390
	g=s(i)
	j=i
390	continue
	if(j.eq.k) go to 450
	s(j)=s(k)
	s(k)=g
	if(nv.eq.0) go to 410
	do 400 i=1,n
	q=v(i,j)
	v(i,j)=v(i,k)
400	v(i,k)=q
410	if(nu.eq.0) go to 430
	do 420 i=1,n
	q=u(i,j)
	u(i,j)=u(i,k)
420	u(i,k)=q
430	if(n.eq.np) go to 450
	do 440 i=n1,np
	q=a(j,i)
	a(j,i)=a(k,i)
        a(k,i)=q
440     continue
450	continue
c
c  back transformation
	if(nu.eq.0) go to 510
	do 500 kk=1,n
	k=n1-kk
	if(b(k).eq.0.0) go to 500
        q=-a(k,k)/cabs(a(k,k))
	do 460 j=1,nu
460	u(k,j)=q*u(k,j)
	do 490 j=1,nu
        q=(0.0,0.0)
	do 470 i=k,m
470     q = q + conjg(a(i,k))*u(i,j)
        q = q/(cabs(a(k,k))*b(k))
	do 480 i=k,m
480	u(i,j)=u(i,j)-q*a(i,k)
490	continue
500	continue
510	if(nv.eq.0) go to 570
	if(n.lt.2) go to 570
	do 560 kk=2,n
	k=n1-kk
	k1=k+1
	if(c(k1).eq.0.0) go to 560
        q=-conjg(a(k,k1))/cabs(a(k,k1))
	do 520 j=1,nv
520	v(k1,j)=q*v(k1,j)
	do 550 j=1,nv
        q=(0.0,0.0)
	do 530 i=k1,n
530	q=q+a(k,i)*v(i,j)
        q=q/(cabs(a(k,k1))*c(k1))
	do 540 i=k1,n
540     v(i,j)=v(i,j)-q*conjg(a(k,i))
550	continue
560	continue
570	return
	end
c******************************************************************************
        subroutine vsb(a,m,n,s,v)
        complex a(m,n),v(n-1,n-1)
c
c  returns result in a(i,1)
c
        real s(n-1)
        do 10 i=1,n-1
   10   a(i,n) = a(i,n)/s(i)
        do 20 j=1,n-1
        a(j,1)=0.0
        do 20 i=1,n-1
   20   a(j,1)=a(j,1)+v(j,i)*a(i,n)
        return
        end


