C----------------------------------------------------------------------
C Copyright (c) 2008, Jens Rademacher
C All rights reserved.
C Copyright details (BSD License) is included at the end of this code
C 
C Modified, with permission, by Matthew Smith, Microsoft Research Ltd.
C for the paper "Explaining spatiotemporal dynamics in lambda-omega
C reaction-diffusion equations using absolute stability spectra of
C wavetrains" 
C by Matthew J. Smith, Jens D.M. Rademacher and Jonathan A. Sherratt.
C SUBMITTED MANUSCRIPT.
C
C----------------------------------------------------------------------

      SUBROUTINE STPNT(NDIM,U,PAR,T) 
C 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      integer k,i
      DIMENSION U(NDIM),PAR(*),V(4)
      DIMENSION UA(NDIM),PARA(13)

	it=1

c---------parameters------
c        par(1) = w0  
c        par(2) = w1           
c        par(3) = R "amplitude"
c        par(4) = real lambda
c        par(5) = imag lambda
c        par(6) = real nu_1 and real nu_2
c        par(7) = imag nu_1
c        par(8) = morse index, but not used here
c        par(9) = eta = (imag nu_1 - imag nu_2)
c        par(10) = real nu_3
c        par(11) = imag nu_3
c        par(12) = real nu_4
c        par(13) = imag nu_4

        if (it.eq.0) then  
c - this case is not used here    
c-----------------
      else
c     This assumes that the triple point is the 
c	first entry in the LO2.dat file
      open(28, FILE='LO2.dat', STATUS='OLD')
      read(28,*)                ! ignore line
      read(28,*) dum, (ua(i),i=1,16)
      do i=1,18
         read(28,*)
      enddo
      read(28,*) (para(i),i=1,13)
      close(28)	

c     this swaps 1 and 4
      u(1)=ua(13)
      u(2)=ua(14)
      u(3)=ua(15)
      u(4)=ua(16)
      u(5)=ua(5)
      u(6)=ua(6)
      u(7)=ua(7)
      u(8)=ua(8)
      u(9)=ua(9)
      u(10)=ua(10)
      u(11)=ua(11)
      u(12)=ua(12)
      u(13)=ua(1)
      u(14)=ua(2)
      u(15)=ua(3)
      u(16)=ua(4)

c     This swaps 2 and 4
c      u(1)=ua(1)
c      u(2)=ua(2)
c      u(3)=ua(3)
c      u(4)=ua(4)
c      u(5)=ua(13)
c      u(6)=ua(14)
c      u(7)=ua(15)
c      u(8)=ua(16)
c      u(9)=ua(9)
c      u(10)=ua(10)
c      u(11)=ua(11)
c      u(12)=ua(12)
c      u(13)=ua(5)
c      u(14)=ua(6)
c      u(15)=ua(7)
c      u(16)=ua(8)

      do i=1,5
         par(i)=para(i)
      enddo

c     This swaps 1 and 4
      par(6)=para(12)
      par(7)=para(13)
      par(8)=para(6)
      par(12)=para(6)
      par(13)=para(7)
      par(10)=para(10)
      par(11)=para(11)
      par(9)=(para(7)+para(9))-para(13)


c     This swaps 2 and 4
c      par(6)=para(6)
c      par(7)=para(7)
c      par(8)=para(6)
c      par(12)=para(6)
c      par(13)=para(7)+para(9)
c      par(10)=para(10)
c      par(11)=para(11)
c      par(9)=-(para(7))+para(13)      

c compute second eigenvector
      v(1) = u(1) - par(9)*u(7)
      v(2) = u(2) - par(9)*u(8)
      v(3) = u(3) + par(9)*u(5)
      v(4) = u(4) + par(9)*u(6)
      do i=1,4
         u(i+4) = v(i)
      enddo
      endif

      RETURN 
      END 

c----------------------------------------
      SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) 

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NDIM),PAR(*),F(NDIM),V(6)

c     Structure:   U = (real1, real2, imaginary1, im2)
c     This program evaluates the matrix Script D. This needs
c     four equations each because of the real and complex parts
c     to the two eigenvectors
c     The equations are different in the case of the second set of 
c     four equations because these correspond to equation (4.8)
c     in Rademacher et al (2008).

      nus = 4                   ! number of nus computed
      it=1
      do i=0,nus-1
         call DISP(V, i, U, PAR,NDIM,IT)
         do j=1,4
            F(j+i*4) = V(j)
         end do
      enddo

c      call myindex(par)         ! compute current index and output a change

      RETURN 
      END 

c----------------------------------------
      subroutine DISP(V, INU, U, PAR, NDIM,it)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      DIMENSION V(6),U(NDIM),PAR(*)

c lambda omega parameters
      w0     = par(1)           !  omega_0
      w1     = PAR(2)           !  omega_1
      R      = PAR(3)           !  amplitude
      FR	 = (1-R**2)     !  'lambda' function
      sqFR	 = FR**(0.5)    !  square root of 'lambda' function
      dFR	 = -2*R		  !  derivative of 'lambda' function
      dGR	 = -2*w1*R		  !  derivative of 'omega' function
	
c spectral parameters
      PLAR = PAR(4)             ! real part of eigenvalue \lambda
      PLAI = PAR(5)             ! imaginary part of eigenvalue
      if(inu.ne.1) then
         PNUR = PAR(6+2*inu)    ! real part of both \nu's
         PNUI = PAR(7+2*inu)    ! imaginary part of first \nu
      else
         pnur = par(6)          ! \nu with same real part as first
         ga   = par(9)
	if (it.eq.0) then
         pnui = par(7)          ! imaginary part of first
	else
	   pnui=par(7)+ga
	endif
      endif

c matrix entries shorthands:
      a1r = (pnur**2-pnui**2)+R*dFR-plar
      a1i = 2*pnui*pnur-plai
      a2r = -2*R*sqFR*pnur
      a2i = -2*R*sqFR*pnui
      a3r = 2*sqFR*pnur/R+dGR
      a3i = 2*sqFR*pnui/R
      a4r = (pnur**2-pnui**2)-plar
      a4i = 2*pnui*pnur-plai

      v1r = U(1 + INU*4)
      v2r = U(2 + INU*4)
      v1i = U(3 + INU*4)
      v2i = U(4 + INU*4)
c real parts:
      V(1) = a1r*v1r - a1i*v1i + a2r*v2r - a2i*v2i
      V(2) = a3r*v1r - a3i*v1i + a4r*v2r - a4i*v2i
c imaginary parts:      
      V(3) = a1i*v1r + a1r*v1i + a2i*v2r + a2r*v2i
      V(4) = a3i*v1r + a3r*v1i + a4i*v2r + a4r*v2i

c desingularized:
      if((inu.eq.1).and.(it.eq.0)) then
c     These first 8 equations correspond to
c     D[(2nu+igamma)+c]
         b1r=2*pnur
         b2r=-2*R*sqFR
         b3r=2*sqFR/R
         b4r=2*pnur

         b1i=2*pnui+ga
         b2i=0
         b3i=0
         b4i=2*pnui+ga
c     These next four correspond to (u+igammav)
         c1r=u(1)-ga*v1i
         c2r=u(2)-ga*v2i
         c1i=u(3)+ga*v1r
         c2i=u(4)+ga*v2rc
c
         V(1)=V(1)+b1r*c1r+b2r*c2r-b1i*c1i-b2i*c2i
         V(2)=V(2)+b3r*c1r+b4r*c2r-b3i*c1i-b4i*c2i
         V(3)=V(3)+b1i*c1r+b2i*c2r+b1r*c1i+b2r*c2i
         V(4)=V(4)+b3i*c1r+b4i*c2r+b3r*c1i+b4r*c2i

      endif      
      RETURN 
      END 

c----------------------------------------
      subroutine myindex(par)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      DIMENSION PAR(*)

      tol = 0.000001

      if(par(6)-par(10).gt.0) then
         if(par(6)-par(12).gt.0) then
            ind = 1
         else
            ind = 2
         endif
      else
         if(par(6)-par(12).gt.0) then
            ind = 2
         else
            if(abs(par(10)-par(12)).gt.tol) then
               ind = 3
            else
               ind = 1
            endif
         endif
      endif
      
c      if(ind.ne.par(8)) then 
         par(8) = ind
c      endif
      return
      end

c----------------------------------------
      SUBROUTINE ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,FI,IJAC,DINT)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NDIM),UOLD(NDIM),UDOT(NDIM),UPOLD(NDIM)
      DIMENSION FI(NINT),ICP(*),PAR(*)
      INTEGER j

c     number of eigenvectors
      NUS=4

c     normalization and phase conditions

      do i=0,NUS-1
c     re(u1), im(u1)
         FI(1+i*2)= -1.0
         FI(2+i*2)=0.0
         do j=1,2
            FI(1+i*2) = FI(1+i*2) + u(j+i*4)*u(j+i*4) + u(j+2+i*4)*u(j+2
     $           +i*4)
            FI(2+i*2) = FI(2+i*2) + UOLD(j+i*4)*U(j+2+i*4) - UOLD(j+2+i
     $           *4)*U(j+i*4)
         end do
      enddo
c other eigenvector
c       FI(3) = 0.0
c       FI(4) = 0.0
c       gam = par(9)
c       do j=5,6
c          VRold = uold(j)
c          VIold = uold(j+2)
c          VR = u(j)
c          VI = u(j+2)
C     im(v1)
c          FI(3) = FI(3) + VIold*u(j-4) - VRold*u(j-2) + VI*uold(j-4) -
c     $         VR*uold(j-2) + GAM*( VRold*VR + VIold *VI )
C     re(v1)
c          FI(4) = FI(4) + VRold*U(j-4) + VIold*u(j-2) - VR*uold(j-4) -
c     $         VI*uold(j-2) - GAM*(VIold*VR - VI*VRold )
c       end do
      RETURN 
      END 

c----------------------------------------
      SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,FB,IJAC,DBC) 
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      DIMENSION PAR(*),ICP(*),U0(NDIM),U1(NDIM),FB(NBC)
      INTEGER j
C     periodic bc
       do j=1,NDIM
          FB(j) = U0(j) - U1(j)
       end do
c
      RETURN 
      END 

      SUBROUTINE FOPT 
      RETURN 
      END 
C 
      SUBROUTINE PVLS(NDIM,PAR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION PAR(*)
      PAR(15)=PAR(6)-PAR(10)
      PAR(16)=PAR(6)-PAR(12)
      RETURN 
      END 

C----------------------------------------------------------------------
C Copyright (c) 2008, Jens Rademacher
C All rights reserved.
C Modified, with permission, by Matthew Smith, Microsoft Research Ltd.
C
C Redistribution and use in source and binary forms, with or without
C  modification, are permitted provided that the following conditions 
C are met:
C 1) Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer. 
C 2) Redistributions in binary form must reproduce the above copyright 
C notice, this list of conditions and the following disclaimer in the 
C documentation and/or other materials provided with the distribution. 
C 3) Neither the name of the contributors, nor the names of their 
C organizations, may be used to endorse or promote products derived 
C from this software without specific prior written permission. 
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.