c
c
       program nonplanar
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
c
c     stability analysis for strong, slow, parallel, ideal MHD shocks 
c
c     Ref: Edelman, M.A. 1989, Astrofizika, 31 #2, 4071
c
c     Version: Janueary 9, 2019
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
c
      parameter (neq=11,i1=50,i2=50)
      dimension y(neq),yprime(neq),yscal(neq)
      common/steady/chib,chic,alpha,beta
      common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
      dimension q(4,2),qf(4,8),qdf(4,4),qnot(4)
      dimension dr(100),di(100),br(100),bi(100)
c
      data eps,max_step/1.0e-08,1000000/
c
c        1         2         3         4         5         6         7 2
c
c     power law cooling function: lambda = C*rho**(beta-alpha)*P**alpha
c
      write(6,*) 'Lambda = constant x rho^beta x temp^alpha'
      write(6,*) ' '

      write(6,*) 'Enter beta'
      read(5,*) beta 

      write(6,*) 'Enter alpha, delta alpha, iter'
      read(5,*)   alpha,delta_alpha,max_alpha
c
      write(6,*) 'Enter xcut--velocity cut-off at bottom of
     &             cooling layer'
      read(5,*)   xcut

      xmin        = - 0.01 * xcut
c
c     chib:  magnetic parameter = 0.5 * [U(mag)/U(flow)]
c
      write(6,*) 'Enter magnetic parameter'
      read(5,*)   chib
c
      if(chib.eq.0.0) then
        nmin = 2
        nmhd = 9
      else
        nmin = 4
        nmhd = 11
      end if
c
c     perturbation wave vector (corrugation inverse wavelength)
c
      write(6,*) 'Enter wave, delta kappa, iter '
      read(5,*)  wave,delta_wave,max_kappa
c
c     eigenvalue guesses
c
      write(6,*) 'Enter range for initial eigenvalues'
      write(6,*) '  fr1,fr2,fi1,fi2'
      write(6,*) '  byr1,byr2,byi1,byi2'
      read(5,*)  q(1,1),q(1,2),q(2,1),q(2,2)
      read(5,*)  q(3,1),q(3,2),q(4,1),q(4,2)
c
c    max_alpha: number of cooling function temperature exponents to
c               consider
c
      do n_alpha = 1 , max_alpha
c  
c     calculate steady-state shock structure
c 
      call equilibrium
c
c     max_kappa: number of corrugation wave vectors to consider
c
      do n_kappa = 1 , max_kappa
c
c ======================================================================
c
c  start iteration loop to find eigenvalues for given cooling function 
c  exponent and corrugation wave vector
c
c ======================================================================
c
 111  continue

      do 200 njacobi = 1 , 2*nmin+1

        go to (1,2,3,4,5,6,7,8,9) njacobi

 1      continue
          fr   = q(1,1)
          fi   = 0.5*(q(2,1)+q(2,2))
          byr1 = 0.5*(q(3,1)+q(3,2))
          byi1 = 0.5*(q(4,1)+q(4,2))
          go to 10
 2      continue
          fr   = q(1,2)
          fi   = 0.5*(q(2,1)+q(2,2))
          byr1 = 0.5*(q(3,1)+q(3,2))
          byi1 = 0.5*(q(4,1)+q(4,2))
          go to 10
 3      continue
          fr   = 0.5*(q(1,1)+q(1,2))
          fi   = q(2,1)
          byr1 = 0.5*(q(3,1)+q(3,2))
          byi1 = 0.5*(q(4,1)+q(4,2))
          go to 10
 4      continue
          fr   = 0.5*(q(1,1)+q(1,2))
          fi   = q(2,2)
          byr1 = 0.5*(q(3,1)+q(3,2))
          byi1 = 0.5*(q(4,1)+q(4,2))
          go to 10
 5      continue
          fr   = 0.5*(q(1,1)+q(1,2))
          fi   = 0.5*(q(2,1)+q(2,2))
          byr1 = q(3,1)
          byi1 = 0.5*(q(4,1)+q(4,2))
          go to 10
 6      continue
          fr   = 0.5*(q(1,1)+q(1,2))
          fi   = 0.5*(q(2,1)+q(2,2))
          byr1 = q(3,2)
          byi1 = 0.5*(q(4,1)+q(4,2))
          go to 10
 7      continue
          fr   = 0.5*(q(1,1)+q(1,2))
          fi   = 0.5*(q(2,1)+q(2,2))
          byr1 = 0.5*(q(3,1)+q(3,2))
          byi1 = q(4,1)
          go to 10
 8      continue
          fr   = 0.5*(q(1,1)+q(1,2))
          fi   = 0.5*(q(2,1)+q(2,2))
          byr1 = 0.5*(q(3,1)+q(3,2))
          byi1 = q(4,2)
          go to 10
 9      continue
          fr   = 0.5*(q(1,1)+q(1,2))
          fi   = 0.5*(q(2,1)+q(2,2))
          byr1 = 0.5*(q(3,1)+q(3,2))
          byi1 = 0.5*(q(4,1)+q(4,2))
          go to 10
 10     continue
c
c ======================================================================
c
c    Start shock structure calculation for given frequency and transverse
c    magnetic field vector, By 
c
        call jump_edelman(x,y,yprime,1)

        htry = - 0.01
        do j = 1 , max_step 
          call derivs(x,y,yprime)
          do k = 1, nmhd
            yscal(k) = abs(y(k))+abs(x*yprime(k))+1.0e-30
          end do
          call rkqc(y,yprime,nmhd,x,htry,eps,yscal,hdid,hnext)
          htry = hnext
          xmax  = xcut - x
          if(htry.lt.xmax) htry = xmax
          if(x.le.xcut) go to 100
        end do
c
c    End of shock structure solution
c
c ======================================================================
c
 100  continue
c
        if (njacobi.eq.2*nmin+1) then
c
          qnot(1) = - y(6)
          qnot(2) = - y(7)

          qnot(3) = - y(10)
          qnot(4) = - y(11)
c
        else
c
          qf(1,njacobi) = y(6)
          qf(2,njacobi) = y(7)

          qf(3,njacobi) = y(10)
          qf(4,njacobi) = y(11)

          nprint=1
          if(nprint.eq.1) then 
            write(6,*) ' '
            write(6,20) alpha,wave,chic,chib
            write(6,20) fr,fi,byr1,byi1,pzr,pzi
            write(6,20) x,y(1),(qf(ii,njacobi),ii=1,4)
 20         format(1x,1p8e10.2)
          end if
c
        end if
c
 200  continue
c
        do i = 1 , nmin
          do j = 1 , nmin
            qdf(i,j) = (qf(i,2*j)-qf(i,2*j-1))/(q(j,2)-q(j,1))
          end do
        end do
c
        call gaussj(qdf,nmin,4,qnot,1)

        on_key = 0
        if(abs(qnot(1)).gt.1.0e-12) on_key = 1
        if(abs(qnot(2)).gt.1.0e-12) on_key = 1

        q_max = 0.0
        do i = 1 , nmin 
          q(i,1) = 0.5 * (q(i,1)+q(i,2))
          if(i.eq.1 .or. i.eq.2) then
            ratio = abs(qnot(i)/q(i,1))
            if(ratio.gt.q_max) q_max = ratio
          end if
        end do

        q_max = amin1( q_max , 1.0e02 )

        if (q_max.lt.0.2) then
          q_max = 1.0
        else
          q_max = 0.2 / q_max
        end if 
        
        do i = 1 , nmin
          qnot(i) = q_max * qnot(i)
          q(i,2) = q(i,1) + 2.0 * qnot(i)
        end do

        if(on_key.eq.1)  go to 111
c
c ======================================================================
c
c  end of iteration loop to find eigenvalues for given cooling function 
c  exponent and corrugation wave vector
c
c ======================================================================
c
        cut=1.0e-20
        if(abs(q(1,1)).le.cut) q(1,1)=cut
        if(abs(q(2,1)).le.cut) q(2,1)=cut
        if(abs(pzr).le.cut) pzr=cut
        write(6,451) alpha,beta,wave,chic,chib
        write(61,451) alpha,beta,wave,chib,
     1   q(1,1),q(2,1),q(3,1),q(4,1)
        write(62,451) alpha,beta,wave,chib,
     1   q(1,1),q(2,1),pzr,pzi
 451    format(1p8e11.3)
c
c ======================================================================
c
c    Calculate and output shock structure for given frequency and By vector 
c
       go to 459

        call jump_edelman(x,y,yprime,1)
        write(63,452) x,(y(i),i=2,11)
 452    format(1p12e10.2)
        htry = - 0.01
        do j = 1 , max_step 
          call derivs(x,y,yprime)
          do k = 1, nmhd
            yscal(k) = abs(y(k))+abs(x*yprime(k))+1.0e-30
          end do
          call rkqc(y,yprime,nmhd,x,htry,eps,yscal,hdid,hnext)
          write(63,452) x,(y(i),i=2,11)
          htry = hnext
          xmax  = xcut - x
          if(htry.lt.xmax) htry = xmax
          if(x.le.xcut) go to 300
        end do

 459  continue
c
c    End of shock structure solution
c
c ======================================================================
c
 300    continue
c
        wave = wave + delta_wave

        q(1,2) = 1.01 * q(1,1)
        q(2,2) = 1.01 * q(2,1)
        q(3,2) = 1.01 * q(3,1)
        q(4,2) = 1.01 * q(4,1)
       
        if ( n_kappa .eq. 1 ) then
          dr( n_alpha ) = q(1,1)
          di( n_alpha ) = q(2,1)
          br( n_alpha ) = q(3,1)
          bi( n_alpha ) = q(4,1)
        end if
 
        end do 

        wave  = wave - max_kappa * delta_wave
        alpha = alpha + delta_alpha
        
        q(1,1) = dr( n_alpha )
        q(1,2) = 1.01 * q(1,1)
        q(2,1) = di(n_alpha) 
        q(2,2) = 1.01 * q(2,1)
        q(3,1) = br(n_alpha) 
        q(3,2) = 1.01 * q(3,1)
        q(4,1) = bi(n_alpha) 
        q(4,2) = 1.01 * q(4,1)

        end do

      end
c=============================================================
      subroutine equilibrium
      parameter (neq=11)
      dimension y(neq),yprime(neq),yout(neq)
      common/steady/chib,chic,alpha,beta

      data max_step/100000/

      chic=1000.0
 100  continue
c
      call jump_edelman(x,y,yprime,0)
c
      do j = 1 , max_step
        call derivs(x,y,yprime)
        xstep=-0.01*abs(y(1)/yprime(1))
        call rk4(y,yprime,1,x,xstep,yout)
        y(1) = yout(1)
        x    = x + xstep 
        if(y(1).ge.-2.5e-08) go to 1000
      end do
c
 1000 continue
c
      frac=(1.0-x)
      if(abs(x).lt.1.0e-12) go to 2000
      chic=chic*frac
      go to 100
c
 2000 continue
c
      return
      end
c====================================================================
      function dv(x1)
      common/steady/chib,chic,alpha,beta
      data alpha1,beta1,ratio/1.0,1.0,1.0/

      cool = chic*(-x1)**(alpha-beta)*(1.0+x1)**alpha
      cool = cool+chic*ratio*(-x1)**(alpha1-beta1)*(1.0+x1)**alpha1
      dv = -2.0*cool/(5.0+8.0*x1)

      return
      end
c====================================================================
      SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT)
      PARAMETER (neq=11)
      DIMENSION Y(neq),DYDX(neq),YOUT(neq),Y1(neq),Y2(neq),Y3(neq)
C
      HH=H/2.0
      H6=H/6.0
      XH=X+HH
C
      DO 11 I=1,N
        YOUT(I)=Y(I)+HH*DYDX(I)
11    CONTINUE
      CALL DERIVS(XH,YOUT,Y1)
      DO 12 I=1,N
        YOUT(I)=Y(I)+HH*Y1(I)
12    CONTINUE
      CALL DERIVS(XH,YOUT,Y2)
      DO 13 I=1,N
        YOUT(I)=Y(I)+H*Y2(I)
13    CONTINUE
      CALL DERIVS(X+H,YOUT,Y3)
C
      DO 14 I=1,N
        YOUT(I)=Y(I)+H6*(DYDX(I)+2.0*(Y1(I)+Y2(I))+Y3(I))
14    CONTINUE
C
      RETURN
      END
c==================================================================
      subroutine jump_edelman(x,y,yprime,n_order)
c
c     n_order -- 0, steady state shock jump conditions
c     n_order -- 1, perturbed shock jump conditions
c
      parameter (neq=11)
      real kzr,kzi
      dimension y(neq),yprime(neq)
      common/steady/chib,chic,alpha,beta
      common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
      complex zimag,zfreq,zkappa
      complex deltab,zdby,zdrho,zdvy
      complex zb1,zb2
c
      zimag = cmplx(0.0,1.0)
      zdby  = cmplx(byr1,byi1)
      zfreq = cmplx(fr,fi)
c
c  zeroth order jump conditions
c
      x    =  1.00
      y(1) = -0.25

      if (n_order.eq.0) then 
c
c  zeroth order jump condition branch
c
        return
c
      else
c
c  perturbed shock jump condition branch
c
c 
c  Require:  kzr < 0 for outward propagation
c            kzi > 0 for decaying wave at infinity
c
        if (chib.le.0.5) then 

          b1r = 0.0
          b2r = 0.0
          b1i = 0.0
          b2i = 0.0
          byr1 = 0.0
          byi1 = 0.0
          pzr = 0.0
          pzi = 0.0

        else 

          call kappa(kzr,kzi)
          write(6,*) kzr,kzi,' quadratic'
cc
          call kappa_polynomial(kzr,kzi)
c
          stop
ccc          call kappa_edelman(kzr,kzi)
cc
          pzr  = kzr
          pzi  = kzi
          zkappa = cmplx(kzr,kzi)

          ckr  = cos(kzr)*exp(-kzi)
          cki  = sin(kzr)*exp(-kzi)

          bkr  =  (kzr*byr1+kzi*byi1)/(kzr**2+kzi**2)
          bki  =  (kzr*byi1-kzi*byr1)/(kzr**2+kzi**2)

          b1r  =  wave*bkr
          b1i  =  wave*bki

          zb1  =  wave*(zdby/zkappa)

          bkr  =  fi*pzr-fr*pzi
          bki  = -fr*pzr-fi*pzi

          b2r  = (byr1*bkr-byi1*bki)/(kzr**2+kzi**2)
          b2i  = (byi1*bkr+byr1*bki)/(kzr**2+kzi**2)
    
          zb2   = zimag*zfreq*(zdby/zkappa)

        end if
c
c
        fin  = 1.0/(fr**2+fi**2)
c
c
c  y(2),y(3)--density perturbation
        y(2) =  - b1r
        y(3) =  - b1i
c  y(4),y(5)--tangential v perturbation (y-component of velocity)
        y10 = (+3.0*wave*fi*fin+4.0*(1.0-2.0*chib)*byr1)/(1.0-8.0*chib)
        y11 = (+3.0*wave*fr*fin+4.0*(1.0-2.0*chib)*byi1)/(1.0-8.0*chib)
        deltab = (+3.0*zimag*wave/zfreq+4.0*(1.0-2.0*chib)*zdby)
     &             /(1.0-8.0*chib)
        y(4) = y10 + 4.0*b2r
        y(5) = y11 + 4.0*b2i
c  y(6),y(7)--perpendicular v perturbation (z-component of velocity)
        y(6) = - 3.0 
        y(7) =   0.0 
c  y(8),y(9)--pressure perturbation
        y(8) = + 2.0 - b1r 
        y(9) =       - b1i
c  y(10),y(11)--tangential B-field perturbation (y-component of B)
        y(10) =  y10 
        y(11) =  y11 
c
c
        return

      end if

      end
c==================================================================
c=======================================================================
      subroutine derivs(x,y,yprime)
c
c
c     f(1): equilibrium momentum equation
c     f(2),f(3): linearized continuity equation
c     f(4),f(5): linearized y-component of momentum equation
c     f(6),f(7): linearized z-component of momentum equation
c     f(8),f(9): linearized energy equation
c     f(10),f(11): linearized y-component of magnetic equation
c
c     note--the z-component of the magnetic equation is algebraic and 
c           is used to eliminate the z-component of the magnetic field 
c           perturbation
c
      parameter (neq=11)
      dimension y(neq),yprime(neq),f(neq)
      common/steady/chib,chic,alpha,beta
      common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
c
      data alpha1,beta1,ratio/1.0,1.0,1.0/
      data gamma,gmin1/1.66666667,0.66666667/
c
c
      frnew = fr
      finew = fi
c
cc      fr =  finew
cc      fi = -frnew
c
c
      v0   = y(1)
      rr   = y(2)
      ri   = y(3)
      yr   = y(4)
      yi   = y(5)
      zr   = y(6)
      zi   = y(7)
      pr   = y(8)
      pi   = y(9)
      byr  = y(10)
      byi  = y(11)

      fin  = 1.0/(fr**2+fi**2)

      bzr  = v0*wave*fin*(fi*(yr-byr)-fr*(yi-byi))
      bzi  = v0*wave*fin*(fr*(yr-byr)+fi*(yi-byi))

      p0   = 1.0+v0
      dvdr = dv(v0)
      dpdr = dvdr
      cool  = gmin1*chic*(p0)**alpha*(-v0)**(alpha-beta)
      cool1 = gmin1*chic*ratio*(p0)**alpha1*(-v0)**(alpha1-beta1)
c
      f(1) =  dvdr
      f(2) = -(rr*fr-ri*fi+x*dvdr/v0)+wave*v0*yi
      f(3) = -(ri*fr+rr*fi)          -wave*v0*yr
      f(4) = -(yr*fr-yi*fi)-dvdr*yr-wave*p0*pi-wave*x*dpdr*fi*fin
     1        +2.0*chib*wave*bzi
      f(5) = -(yr*fi+yi*fr)-dvdr*yi+wave*p0*pr-wave*x*dpdr*fr*fin
     1        -2.0*chib*wave*bzr
      f(6) = -(zr*fr-zi*fi-x*dvdr/v0)-dvdr*rr+dpdr*pr
     1        -2.0*dvdr*zr
      f(7) = -(zr*fi+zi*fr)-dvdr*ri+dpdr*pi
     1        -2.0*dvdr*zi
      f(8) = -p0*(pr*fr-pi*fi-gamma*(rr*fr-ri*fi))
     1        -cool*(x/v0+(alpha-1.0)*pr+(beta-alpha)*rr+fr*fin-zr)
     1        -cool1*(x/v0+(alpha1-1.0)*pr+(beta1-alpha1)*rr+fr*fin-zr)
      f(9) = -p0*(pr*fi+pi*fr-gamma*(rr*fi+ri*fr))
     1        -cool*((alpha-1.0)*pi+(beta-alpha)*ri-fi*fin-zi)
     1        -cool1*((alpha1-1.0)*pi+(beta1-alpha1)*ri-fi*fin-zi)
      f(10) = -dvdr/v0*(yr-byr)+(byr*fr-byi*fi)/v0
      f(11) = -dvdr/v0*(yi-byi)+(byi*fr+byr*fi)/v0
c
      yprime(1)  =  f(1)

      yprime(2)  = (f(2)-f(6)-f(8)/v0)/(v0+gamma*p0)
      yprime(3)  = (f(3)-f(7)-f(9)/v0)/(v0+gamma*p0)

      yprime(4)  = (f(4)-2.0*chib*f(10))/(v0-2.0*chib)
      yprime(5)  = (f(5)-2.0*chib*f(11))/(v0-2.0*chib)

      yprime(6)  =  f(2)/v0-yprime(2)
      yprime(7)  =  f(3)/v0-yprime(3)

      yprime(8)  =  f(8)/(v0*p0)+gamma*yprime(2)
      yprime(9)  =  f(9)/(v0*p0)+gamma*yprime(3)

      yprime(10) =  yprime(4)-f(10)
      yprime(11) =  yprime(5)-f(11)
c
c
      fr =  frnew
      fi =  finew
c
c
      return
      end
c=====================================================================
      SUBROUTINE GAUSSJ(A,N,NP,B,M)
      PARAMETER (NMAX=50)
      DIMENSION A(4,4),B(4),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
      DO 11 J=1,N
        IPIV(J)=0
11    CONTINUE
      DO 22 I=1,N
        BIG=0.
        DO 13 J=1,N
          IF(IPIV(J).NE.1)THEN
            DO 12 K=1,N
              IF (IPIV(K).EQ.0) THEN
                IF (ABS(A(J,K)).GE.BIG)THEN
                  BIG=ABS(A(J,K))
                  IROW=J
                  ICOL=K
                ENDIF
              ELSE IF (IPIV(K).GT.1) THEN
                PAUSE 'Singular matrix'
              ENDIF
12          CONTINUE
          ENDIF
13      CONTINUE
        IPIV(ICOL)=IPIV(ICOL)+1
        IF (IROW.NE.ICOL) THEN
          DO 14 L=1,N
            DUM=A(IROW,L)
            A(IROW,L)=A(ICOL,L)
            A(ICOL,L)=DUM
14        CONTINUE
          DO 15 L=1,M
            DUM=B(IROW)
            B(IROW)=B(ICOL)
            B(ICOL)=DUM
15        CONTINUE
        ENDIF
        INDXR(I)=IROW
        INDXC(I)=ICOL
        IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.'
        PIVINV=1./A(ICOL,ICOL)
        A(ICOL,ICOL)=1.
        DO 16 L=1,N
          A(ICOL,L)=A(ICOL,L)*PIVINV
16      CONTINUE
        DO 17 L=1,M
          B(ICOL)=B(ICOL)*PIVINV
17      CONTINUE
        DO 21 LL=1,N
          IF(LL.NE.ICOL)THEN
            DUM=A(LL,ICOL)
            A(LL,ICOL)=0.
            DO 18 L=1,N
              A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18          CONTINUE
            DO 19 L=1,M
              B(LL)=B(LL)-B(ICOL)*DUM
19          CONTINUE
          ENDIF
21      CONTINUE
22    CONTINUE
      DO 24 L=N,1,-1
        IF(INDXR(L).NE.INDXC(L))THEN
          DO 23 K=1,N
            DUM=A(K,INDXR(L))
            A(K,INDXR(L))=A(K,INDXC(L))
            A(K,INDXC(L))=DUM
23        CONTINUE
        ENDIF
24    CONTINUE
      RETURN
      END
c=======================================================================
      SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT)
      PARAMETER (FCOR=.0666666667,
     *    ONE=1.,SAFETY=0.9,ERRCON=6.E-4,neq=11)
      DIMENSION Y(neq),DYDX(neq),YSCAL(neq),YTEMP(neq),YSAV(neq)
     1 ,DYSAV(neq)

      PGROW=-0.20
      PSHRNK=-0.25
      XSAV=X
      DO 11 I=1,N
        YSAV(I)=Y(I)
        DYSAV(I)=DYDX(I)
11    CONTINUE
      H=HTRY
1     HH=0.5*H
      CALL RK4(YSAV,DYSAV,N,XSAV,HH,YTEMP)
      X=XSAV+HH
      CALL DERIVS(X,YTEMP,DYDX)
      CALL RK4(YTEMP,DYDX,N,X,HH,Y)
      X=XSAV+H
      IF(X.EQ.XSAV)PAUSE 'Stepsize not significant in RKQC.'
      CALL RK4(YSAV,DYSAV,N,XSAV,H,YTEMP)
      ERRMAX=0.
      DO 12 I=1,N
        YTEMP(I)=Y(I)-YTEMP(I)
        ERRMAX=MAX(ERRMAX,ABS(YTEMP(I)/YSCAL(I)))
12    CONTINUE
      ERRMAX=ERRMAX/EPS
      IF(ERRMAX.GT.ONE) THEN
        H=SAFETY*H*(ERRMAX**PSHRNK)
        GOTO 1
      ELSE
        HDID=H
        IF(ERRMAX.GT.ERRCON)THEN
          HNEXT=SAFETY*H*(ERRMAX**PGROW)
        ELSE
          HNEXT=4.*H
        ENDIF
      ENDIF
      DO 13 I=1,N
        Y(I)=Y(I)+YTEMP(I)*FCOR
13    CONTINUE
      RETURN
      END
c========================================================================
      subroutine kappa_edelman (kzr,kzi)
      real kzr,kzi
      common/steady/chib,chic,alpha,beta
      common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
c
      vpre = -1.0

      wr =   2.0 * fi / ( wave * vpre )
      wi = - 2.0 * fr / ( wave * vpre )

      w1 =   wr * wr - wi * wi
      w2 =   wr * wi

      r1 = ( 2.0 * chib - 1.0 ) 
      r2 = ( 2.0 * chib ) / ( 2.0 * chib - 1.0 ) * * 2

      wrchi =  0.5 * wr / r1 
      wichi =  0.5 * wi / r1 

      b1 = ( 0.25 * w1 - r1 )
      radical = sqrt ( b1 * b1 + 0.25 * w2 * w2 )
c
c  real part of kappa
c
      xlp = sqrt ( 0.5 * r2 * ( b1 + radical ) )

      xrpm = 0.5 * ( wr / r1 ) - xlp
      xrpp = 0.5 * ( wr / r1 ) + xlp
c
c  imaginary part of kappa
c
      xipm = 0.5 * ( wi / r1 ) + r2 * 0.25 * w2 / ( xrpm - wrchi )
      xipp = 0.5 * ( wi / r1 ) + r2 * 0.25 * w2 / ( xrpp - wrchi )
c
c
      xrpm = wave * xrpm
      xrpp = wave * xrpp
c
c 
      xipm = wave * xipm
      xipp = wave * xipp
c
c
      write(6,90) fr,fi
      write(6,110) xrpm,xipm,xrpp,xipp
  90  format(1p2e12.4)
 100  format(3h"-",1p2e12.3,4x,1p2e12.4)
 110  format(3h"+",1p2e12.3,4x,1p2e12.4)
c
c  check for outgoing and damping waves at + infinity
c
      if(xrpm.ge.0.0) then
        krpm = -1
      else
        krpm = 1 
      end if
      if(xrpp.ge.0.0) then
        krpp = -1
      else
        krpp = 1 
      end if
c
c
      if(xipm.ge.0.0) then
        kipm = 1
      else
        kipm = -1 
      end if
      if(xipp.ge.0.0) then
        kipp = 1
      else
        kipp = -1 
      end if
c
c
      if (krpm.eq.1 .and. kipm.eq.1) then
        kzr=xrpm
        kzi=xipm
      end if

      if (krpp.eq.1 .and. kipp.eq.1) then
        kzr=xrpp
        kzi=xipp
      end if
c
      return
      end
c========================================================================
      subroutine kappa (kzr,kzi)
      real kzr,kzi,kzrp,kzip,kzrm,kzim 
      common/steady/chib,chic,alpha,beta
      common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
c  
      complex zfreq,zkappa
c
      zfreq = cmplx(fr,fi)
c
      r1 = ( 2.0 * chib - 1.0 ) * wave * * 2 + ( fr * * 2 - fi * * 2 )
      r2 =   2.0 * fr * fi 

      phase = atan(r2/r1)
      if ( r1 .le. 0.0 ) phase = phase + acos( - 1.0 )

      ai =   ( r1 * * 2 + r2 * * 2 ) * * 0.25 * cos( 0.5 * phase )
      ar = - ( r1 * * 2 + r2 * * 2 ) * * 0.25 * sin( 0.5 * phase )

      ar = sqrt( 2.0 * chib ) / ( 2.0 * chib - 1.0 ) * ar
      ai = sqrt( 2.0 * chib ) / ( 2.0 * chib - 1.0 ) * ai

      br = - fi / ( 2.0 * chib - 1.0 ) 
      bi =   fr / ( 2.0 * chib - 1.0 )

      kzrm = br - ar
      kzim = bi - ai
      kzrp = br + ar
      kzip = bi + ai

      k_success = 0

      if (kzrm.le.0.0) then
        if (kzim.ge.0.0) then
          write(6,*) 'good root for - sign'
          kzr = kzrm
          kzi = kzim
          k_success = 1
        end if
      end if

      if (kzrp.le.0.0) then
        if (kzip.ge.0.0) then
          write(6,*) 'good root for + sign'
          kzr = kzrp
          kzi = kzip
          k_success = 1
        end if
      end if

      if (k_success.eq.0) then 
        write(6,*) 'no good root for either sign'
        kzr   = 1.0
        kzi   = 1.0
      end if

      write(6,*) kzr,kzi
      zkappa = cmplx(kzr,kzi)

      return
      end

c=======================================================================

      subroutine kappa_polynomial(kzr,kzi)
c
c     n_order -- 0, steady state shock jump conditions
c     n_order -- 1, perturbed shock jump conditions
c
      parameter (neq=11)
      real kzr,kzi
      dimension y(neq),yprime(neq)
      common/steady/chib,chic,alpha,beta
      common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
      parameter(mroots=4,mroots1=mroots+1,maxm=101)
      COMPLEX A(mroots1),ROOTS(mroots),AD(MAXM)
c
        xmachsq=400.0e00
        velocity=-1.0

        gamma=(5.0/3.0)

        fedr=fr
        fedi=fi
        fr= fedi
        fi=-fedr
        gmach=gamma/xmachsq 

        bsq=(wave*velocity)**2
        fsqr  = fr**2-fi**2
        fsqi  =+2.0*fr*fi
        fcuber= fr*fsqr-fi*fsqi
        fcubei= fr*fsqi+fi*fsqr
        fquadr= fsqr**2-fsqi**2
        fquadi=+2.0*fsqr*fsqi
c
c
c  Dispersion Relation:  fourth-order polynomial
c
c
        c4r=(1.0-2.0*chib)*(1.0-gmach)
        c4i=0.0

        c3r=(2.0-2.0*chib-gmach)*2.0*fr
        c3i=(2.0-2.0*chib-gmach)*2.0*fi

        c2r=-bsq*2.0*chib-bsq*(1.0-2.0*chib)*gmach
     &     +(6.0-2.0*chib-gmach)*fsqr
        c2i=(6.0-2.0*chib-gmach)*fsqi

        c1r=(-2.0*bsq*(2.0*chib+gmach))*fr+4.0*fcuber
        c1i=(-2.0*bsq*(2.0*chib+gmach))*fi+4.0*fcubei

        c0r=fquadr-bsq*(2.0*chib+gmach)*fsqr
        c0i=fquadi-bsq*(2.0*chib+gmach)*fsqi
c
c
        a(5)=cmplx(c4r,c4i)
        a(4)=cmplx(c3r,c3i)
        a(3)=cmplx(c2r,c2i)
        a(2)=cmplx(c1r,c1i)
        a(1)=cmplx(c0r,c0i)
c
        write(6,991) fr,fi,wave,chib
 991    format(1p4e12.4)
c
c
        call ZROOTS(A,mroots,ROOTS,POLISH)
c
        print *,' '
        do i=1,mroots
          write(6,*) roots(i)/velocity
        end do
        print *,' '
c
c
c  Mach number set to infinity
c
c
        c4r=0.0
        c4i=0.0

        c3r=0.0
        c3i=0.0

        c2r=(1.0-2.0*chib)
        c2i=0.0

        c1r=2.0*fr
        c1i=2.0*fi

        c0r=fsqr-2.0*bsq*chib
        c0i=fsqi
c
c
        a(5)=(0.0,0.0)
        a(4)=(0.0,0.0)
        a(3)=cmplx(c2r,c2i)
        a(2)=cmplx(c1r,c1i)
        a(1)=cmplx(c0r,c0i)
c
        write(6,991) ,fr,fi,wave,chib
c
c
        call ZROOTS(A,2,ROOTS,POLISH)
c
        print *,' '
        do i=1,2
          write(6,*) roots(i)/velocity
        end do
        print *,' '
c
c
        fr=fedr
        fi=fedi

      return
      end

      SUBROUTINE ZROOTS(A,M,ROOTS,POLISH)
      PARAMETER (EPS=1.E-8,MAXM=101)
      parameter(mroots=4,mroots1=mroots+1)
      COMPLEX A(mroots1),ROOTS(mroots),AD(MAXM),X,B,C
      LOGICAL POLISH

      DO 11 J=1,M+1
        AD(J)=A(J)
11    CONTINUE
      DO 13 J=M,1,-1
        X=CMPLX(0.0,0.0)
        CALL LAGUER(AD,J,X,EPS,.FALSE.)
        IF(ABS(AIMAG(X)).LE.2.0*EPS**2*ABS(REAL(X))) X=CMPLX(REAL(X),0.)
        ROOTS(J)=X
        B=AD(J+1)
        DO 12 JJ=J,1,-1
          C=AD(JJ)
          AD(JJ)=B
          B=X*B+C
12      CONTINUE
13    CONTINUE
      IF (POLISH) THEN
        DO 14 J=1,M
          CALL LAGUER(A,M,ROOTS(J),EPS,.TRUE.)
14      CONTINUE
      ENDIF
      DO 16 J=2,M
        X=ROOTS(J)
        DO 15 I=J-1,1,-1
          IF(REAL(ROOTS(I)) .LE. REAL(X))GO TO 10
          ROOTS(I+1)=ROOTS(I)
15      CONTINUE
        I=0
10      ROOTS(I+1)=X
16    CONTINUE
      RETURN
      END

c=======================================================================

      SUBROUTINE LAGUER(A,M,X,EPS,POLISH)
      parameter(mroots=4,mroots1=mroots+1,maxm=101)
      COMPLEX A(mroots1),X,DX,X1,B,D,F,G,H,SQ,GP,GM,G2,ZERO
      LOGICAL POLISH
      PARAMETER (ZERO=(0.0,0.0),EPSS=6.E-8,MAXIT=100)

      DXOLD=CABS(X)
      DO 12 ITER=1,MAXIT
        B=A(M+1)
        ERR=CABS(B)
        D=ZERO
        F=ZERO
        ABX=CABS(X)
        DO 11 J=M,1,-1
          F=X*F+D
          D=X*D+B
          B=X*B+A(J)
          ERR=CABS(B)+ABX*ERR
11      CONTINUE
        ERR=EPSS*ERR
        IF(CABS(B) .LE. ERR) THEN
          DX=ZERO
          RETURN
        ELSE
          G=D/B
          G2=G*G
          H=G2-2.0*F/B
          SQ=CSQRT((M-1)*(M*H-G2))
          GP=G+SQ
          GM=G-SQ
          IF(CABS(GP) .LT. CABS(GM)) GP=GM
          DX=M/GP
        ENDIF
        X1=X-DX
        IF(X.EQ.X1)RETURN
        X=X1
        CDX=CABS(DX)
        IF(ITER.GT.6 .AND. CDX.GE.DXOLD) RETURN
        DXOLD=CDX
        IF(.NOT.POLISH)THEN
          IF(CABS(DX) .LE. EPS*CABS(X)) RETURN
        ENDIF
12    CONTINUE
      PAUSE 'too many iterations'
      RETURN
      END

c=======================================================================