      subroutine euler(edim,mdim,iwin,err,ipen,n,x,z,a)
      common /work1/xc(512)
      common /work2/dx(512),dz(512)
      dimension a(n),x(n),z(n)
c      dimension b(7,4),bt(4,7),btb(4,4),binvbt(4,7),e(4),f(4),c(7)
      dimension b(80),bt(80),btb(16),binvbt(80),f(4)
      complex xc
      character*1 sym(7)
      data sym/'c','s','s','p','p','d','d'/
      size=.08
      call spline(a,dx,n,x(2)-x(1))
      call setfft(0,lc,n,dx)
      k=lc/2
      do 30 j=2,k
      xc(j)=xc(j)+xc(j)
   30 xc(lc+2-j)=0.
      call fork(lc,xc,-1.)
      do 31 j=1,n
   31 dz(j)=-aimag(xc(j))
c      call dplot1(xmin,xmax,zmin,zmax,1.)
c      open(20,file='eov.out',status='unknown',form='formatted')
c      do 60 kk=1,30
cc      print *,kk,mdim,edim
c      if(mdim.eq.3) then
c        edim=kk*0.1
c      else
c        if(kk.eq.2) return
c      endif
      if(err.eq.0.) then
        std=1.e38
      else
        std=100./err
      endif
      if(mdim.eq.4) 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 60 k=1,n-iwin+1
      call esolve(d,k,edim,mdim,iwin,b,bt,btb,binvbt,f,n,a,dx,dz,x)
      if(d.eq.0.0) go to 61
      std=0.0
      do 42 i=1,iwin
      e=(f(1)-x(k+i-1))*dx(k+i-1)+f(2)*dz(k+i-1)+f(3)
      if(mdim.eq.4) then
c        e=a(k+i-1)-e/f(4)
        e=e-f(4)*a(k+i-1)
      else
c        e=a(k+i-1)-e/edim
        e=e-edim*a(k+i-1)
      endif
   42 std=std+e*e
      e=sqrt(std/float(iwin-mdim))
   61 continue
c      write(20,200) x(k+3),a(k+3),dx(k+3),dz(k+3),f
  200 format(5f15.6)
      if(f(2).lt.0.0) go to 60
      if(mdim.eq.4) then
        if(f(4).lt.0.0) go to 60
        if(f(4).gt.3.0) go to 60
      endif
      if(mdim.eq.4) then
cc        std=1./(f(2)**2*f(4)*std)
cc        std=1./(f(4)*std)
         if(btb(6).lt.0.) go to 60
         e1=e*sqrt(btb(6))
         edim=f(4)
cc         std=f(2)/(f(4)*sqrt(btb(6)))
c         std=1./(f(4)*e)
c         write(20,200) f(4),f(1),f(2),e
      else
cc        std=1./(f(2)**2*edim*std)
cc        std=1./(edim*std)
         if(btb(5).lt.0.) go to 60
         e1=e*sqrt(btb(5))
cc         std=f(2)/(edim*sqrt(btb(5)))
c         std=1./(edim*e)
c         write(20,200) edim,f(1),f(2),e
      endif
      std=f(2)/e1
c      write(20,200) edim,f(1),f(2),e1,std,e
      if(std.lt.err) go to 60
cc      if(std.gt.err) go to 60
c      std=f(2)*(20.+1./e)
c      if(std.gt.err+5.) go to 60
c      if(std.lt.err-5.) go to 60
      do 43 i=2,n
      if(f(1).gt.x(i)) go to 43
      go to 44
   43 continue
      go to 60
   44 f(2)=f(2)+z(i-1)+(f(1)-x(i-1))*(z(i)-z(i-1))/(x(i)-x(i-1))
      std=100./std
      if(mdim.eq.4) then
        i=ifix(2.*f(4))+1
        write(20,200) f(1),f(2),f(4),std
      else
        i=ifix(2.*edim)+1
        write(20,200) f(1),f(2),edim,std
      endif
      call vchar(f(1),f(2),sym(i),1,ipen,size,0.,0.,0.)
c      print *, f(1), f(2), f(3), std
   60 continue
      return
      end
c
      subroutine esolve(d,k,edim,md,iw,b,bt,btb,binvbt,f,n,a,dx,dz,x)
      dimension b(iw,md),bt(md,iw),btb(md,md)
      dimension binvbt(md,iw),f(md),c(20)
      dimension a(n),dx(n),dz(n),x(n)
      dimension l(4),m(4)
      do 40 i=1,iw
      b(i,1)=dx(k+i-1)
      b(i,2)=dz(k+i-1)
      b(i,3)=1.0
      if(md.eq.4) b(i,4)=-a(k+i-1)
      c(i)=x(k+i-1)*b(i,1)
   40 if(md.eq.3) c(i)=c(i)+edim*a(k+i-1)
      do 41 i=1,iw
      do 41 j=1,md
   41 bt(j,i)=b(i,j)
      call mmult(bt,b,btb,md,iw,md,md,iw,iw,md,md,md)
      call minv(btb,md,d,l,m)
      if(d.eq.0.0) return
      call mmult(btb,bt,binvbt,md,md,iw,md,md,md,iw,md,iw)
      call mmult(binvbt,c,f,md,iw,1,md,iw,iw,1,md,1)
      return
      end
c
C     SUB-ROUTINE TO FIND INVERSE OF MATRIX A. THE INPUT MATRIX A
C     IS DESTROYED DURING COMPUTATION AND REPLACED BY THE
C     RESULTANT INVERSE. SEE IBM'S SSP.
C     A = INPUT MATRIX
C     N = ORDER OF MATRIX A
C     D = RESULTANT DETERMINANT
C     L,M = WORK VECTORS OF LENGTH N
C=====================================================================
C
      SUBROUTINE MINV(A,N,D,L,M)
c      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(1),L(1),M(1)
C        SEARCH FOR LARGEST ELEMENT
      D=1.0
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
   10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20
   15 BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE
C        INTERCHANGE ROWS
      J=L(K)
      IF(J-K) 35,35,25
   25 KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
   30 A(JI)=HOLD
C     INTERCHANGE COLUMNS
   35 I=M(K)
      IF(I-K) 45,45,38
   38 JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
   40 A(JI)=HOLD
C        DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS
C        CONTAINED IN BIGA
   45 IF(BIGA) 48,46,48
   46 D=0.0
      RETURN
   48 DO 55 I=1,N
      IF(I-K) 50,55,50
   50 IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE
C        REDUCE MATRIX
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I-K) 60,65,60
   60 IF(J-K) 62,65,62
   62 KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE
C        DIVIDE ROW BY PIVOT
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J-K) 70,75,70
   70 A(KJ)=A(KJ)/BIGA
   75 CONTINUE
C        PRODUCT OF PIVOTS
      D=D*BIGA
C        REPLACE PIVOT BY RECIPROCAL
      A(KK)=1.0/BIGA
   80 CONTINUE
C        FINAL ROW AND COLUMN INTERCHANGE
      K=N
  100 K=(K-1)
      IF(K) 150,150,105
  105 I=L(K)
      IF(I-K) 120,120,108
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
  110 A(JI)=HOLD
  120 J=M(K)
      IF(J-K) 100,100,125
  125 KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
  130 A(JI)=HOLD
      GO TO 100
  150 RETURN
      END
C======================================================================
C     SUBROUTINE FOR MULTIPLICATION OF MATRIX A BY MATRIX B
C     TO GIVE MATRIX C. A IS L ROWS BY N COLUMNS.
C     B IS N ROWS BY M COLUMNS ,AND C WILL BE L ROWS BY M COLUMNS.
C======================================================================
C
      SUBROUTINE MMULT(A,B,C,L,N,M,NA,MA,NB,MB,NC,MC)
c      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(NA,MA),B(NB,MB),C(NC,MC)
      DO 100 I=1,L
      DO 101 J=1,M
      C(I,J)=0.0
      DO 102 K=1,N
      C(I,J)=C(I,J)+A(I,K)*B(K,J)
  102 CONTINUE
  101 CONTINUE
  100 CONTINUE
      RETURN
      END

