c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
c     File:        mopagrat.f (Fortran 77 source for the mopagrat routine)
c     Author:      Fredrik Jonsson <fj@phys.soton.ac.uk>
c     Date:        December 26, 2005
c     Last change: January 3, 2006
c
c     Description: The mopagrat routine calculates magneto-optically
c                  parametrically amplified field envelopes of forward
c                  and backward travelling wave components inside a Bragg
c                  grating of sinusoidal spatial distribution of the
c                  refractive index.
c
c     Required external routines:
c
c        ccubsolv - For solving the characteristic complex-valued cubic
c                   polynomial equation that appear for the eigenvalues
c                   of the system determinant of the wave propagation.
c
c        zgaussj  - For solving complex-valued linear systems of equations
c                   by Gauss-Jordan elimination.
c
c     Copyright (c) 2006, Fredrik Jonsson
c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
      subroutine mopagrat(alpha,beta,gamma,delta,eta,kappa,ll,z,
     +     a1iplus,a2iplus,a1iminus,a2iminus,
     +     a1famplplus,a1bamplplus,a2famplplus,
     +     a1famplminus,a1bamplminus,a2famplminus,
     +     a1fplus,a1bplus,a2fplus,a1fminus,a1bminus,a2fminus)

      real*8 alpha,beta,gamma,delta,eta,kappa,ll,z
      complex*8 a1iplus,a2iplus,a1iminus,a2iminus
      real*8 a1famplplus,a1bamplplus,a2famplplus
      real*8 a1famplminus,a1bamplminus,a2famplminus
      complex*8 a1fplus,a1bplus,a2fplus,a1fminus,a1bminus,a2fminus,c(3)
      complex*8 bbplus(3,3),xiplus(3,3),lambdaplus(3),ccplus(3)
      complex*8 bbminus(3,3),ximinus(3,3),lambdaminus(3),ccminus(3)
      real*8 tmp
      integer*4 k
      external ccubsolv,zgaussj

c     *********************************************************************
c     Calculate the eigenvalues of the system matrix A of the magneto-
c     optical parametric interaction. This is done by solving the cubic
c     polynomial equation det(A-I*lambda)=0, where I=diag(1,...,1) denotes
c     the identity matrix, or equivalently of the form
c
c             lambda^3+c[2]*lambda^2+c[1]*lambda+c[0]=0,
c
c     by using the external ccubsolv() routine. Here the coefficients are
c     as follows:
c
c             c(1)=c[0],  c(2)=c[1],  c(3)=c[2].
c
c     After the call, the three complex-valued eigenvalues are returned
c     in the vector lambda(1..3).
c     *********************************************************************
      c(1)=cmplx(0.0,((abs(kappa)**2)-((alpha+delta)**2))
     +     *(alpha+delta+beta+gamma)-((eta)**2)*(alpha+delta))
      c(2)=cmplx(((alpha+delta)**2)-((kappa)**2)-((eta)**2),0.0)
      c(3)=cmplx(0.0,-(alpha+delta+beta+gamma))
      call ccubsolv(c,lambdaplus)
      c(1)=cmplx(0.0,((abs(kappa)**2)-((alpha-delta)**2))
     +     *(alpha-delta+beta-gamma)-((eta)**2)*(alpha-delta))
      c(2)=cmplx(((alpha-delta)**2)-((kappa)**2)-((eta)**2),0.0)
      c(3)=cmplx(0.0,-(alpha-delta+beta-gamma))
      call ccubsolv(c,lambdaminus)

c     *********************************************************************
c     Having calculated the eigenvalues of the system matrix A, now
c     proceed with calculating the corresponding eigenvectors. The eigen-
c     vectors are stored in the matrix xi, in which the (j,k):th element
c     is the j:th element of the k:th eigenvector.
c     *********************************************************************
      do 10,k=1,3
         xiplus(1,k)=cmplx(1.0,0.0)
         xiplus(2,k)=cmplx(kappa,0.0)
         xiplus(2,k)=xiplus(2,k)/(cmplx(0.0,alpha+delta)
     +        +lambdaplus(k))
         xiplus(3,k)=(cmplx(0.0,alpha+delta)-lambdaplus(k))
         xiplus(3,k)=xiplus(3,k)+cmplx((kappa**2),0.0)
     +        /(cmplx(0.0,alpha+delta)+lambdaplus(k))
         xiplus(3,k)=xiplus(3,k)*cmplx(0.0,1.0/eta)
         ximinus(1,k)=cmplx(1.0,0.0)
         ximinus(2,k)=cmplx(kappa,0.0)
         ximinus(2,k)=ximinus(2,k)/(cmplx(0.0,alpha-delta)
     +        +lambdaminus(k))
         ximinus(3,k)=(cmplx(0.0,alpha-delta)-lambdaminus(k))
         ximinus(3,k)=ximinus(3,k)+cmplx((kappa**2),0.0)
     +        /(cmplx(0.0,alpha-delta)+lambdaminus(k))
         ximinus(3,k)=ximinus(3,k)*cmplx(0.0,1.0/eta)
 10   continue

c     *********************************************************************
c     Calculate the inverse of the matrix bb by using the Gauss--Jordan
c     elimination routine zgaussj(), as customized from the Numerical
c     Recipes' real-valued version of the same algorithm.
c     *********************************************************************
      tmp=max(realpart(lambdaplus(1)),
     +     realpart(lambdaplus(2)),realpart(lambdaplus(3)))
      do 20,k=1,3
         bbplus(1,k)=xiplus(1,k)
         bbplus(2,k)=xiplus(2,k)*
     +        cexp((lambdaplus(k)-cmplx(tmp,0.0))*cmplx(ll,0.0))
         bbplus(3,k)=xiplus(3,k)
 20   continue
      ccplus(1)=a1iplus
      ccplus(2)=cmplx(0.0,0.0)
      ccplus(3)=a2iplus
      call zgaussj(bbplus,3,3,ccplus,1,1)
      tmp=max(realpart(lambdaminus(1)),
     +     realpart(lambdaminus(2)),realpart(lambdaminus(3)))
      do 30,k=1,3
         bbminus(1,k)=ximinus(1,k)
         bbminus(2,k)=ximinus(2,k)*
     +        cexp((lambdaminus(k)-cmplx(tmp,0.0))*cmplx(ll,0.0))
         bbminus(3,k)=ximinus(3,k)
 30   continue
      ccminus(1)=a1iminus
      ccminus(2)=cmplx(0.0,0.0)
      ccminus(3)=a2iminus
      call zgaussj(bbminus,3,3,ccminus,1,1)

c     *********************************************************************
c     Calculate the logarithmic amplification of the transmitted and
c     reflected signal and idler, in measures of dB.
c     *********************************************************************
      a1fplus=cmplx(0.0,0.0)
      a1bplus=cmplx(0.0,0.0)
      a2fplus=cmplx(0.0,0.0)
      a1fminus=cmplx(0.0,0.0)
      a1bminus=cmplx(0.0,0.0)
      a2fminus=cmplx(0.0,0.0)
      do 40,k=1,3
         tmp=exp(realpart(lambdaplus(k))*ll)
         a1fplus=a1fplus+ccplus(k)*xiplus(1,k)
     +        *cmplx(tmp*cos(imagpart(lambdaplus(k))*ll),
     +        tmp*sin(imagpart(lambdaplus(k))*ll))
         a1bplus=a1bplus+ccplus(k)*xiplus(2,k)
         a2fplus=a2fplus+ccplus(k)*xiplus(3,k)
     +        *cmplx(tmp*cos(imagpart(lambdaplus(k))*ll),
     +        tmp*sin(imagpart(lambdaplus(k))*ll))
         tmp=exp(realpart(lambdaminus(k))*ll)
         a1fminus=a1fminus+ccminus(k)*ximinus(1,k)
     +        *cmplx(tmp*cos(imagpart(lambdaminus(k))*ll),
     +        tmp*sin(imagpart(lambdaminus(k))*ll))
         a1bminus=a1bminus+ccminus(k)*ximinus(2,k)
         a2fminus=a2fminus+ccminus(k)*ximinus(3,k)
     +        *cmplx(tmp*cos(imagpart(lambdaminus(k))*ll),
     +        tmp*sin(imagpart(lambdaminus(k))*ll))
 40   continue
      a1famplplus=cabs(a1fplus/a1iplus)**2
      a1bamplplus=cabs(a1bplus/a1iplus)**2
      a2famplplus=cabs(a2fplus/a1iplus)**2
      a1famplminus=cabs(a1fminus/a1iminus)**2
      a1bamplminus=cabs(a1bminus/a1iminus)**2
      a2famplminus=cabs(a2fminus/a1iminus)**2

c     *********************************************************************
c     Finally calculate the envelopes of forward and backward traveling
c     components of the signal (a1f,a1b) and idler (a2f) fields at spatial
c     location given by z (with 0 < z < L).
c     *********************************************************************
      if (z.ne.ll) then
         a1fplus=cmplx(0.0d0,0.0d0)
         a1bplus=cmplx(0.0d0,0.0d0)
         a2fplus=cmplx(0.0d0,0.0d0)
         a1fminus=cmplx(0.0d0,0.0d0)
         a1bminus=cmplx(0.0d0,0.0d0)
         a2fminus=cmplx(0.0d0,0.0d0)
         do 50,k=1,3
            a1fplus=a1fplus+ccplus(k)*xiplus(1,k)
     +           *cexp(lambdaplus(k)*cmplx(z,0.0))
            a1bplus=a1bplus+ccplus(k)*xiplus(2,k)
     +           *cexp(lambdaplus(k)*cmplx(z,0.0))
            a2fplus=a2fplus+ccplus(k)*xiplus(3,k)
     +           *cexp(lambdaplus(k)*cmplx(z,0.0))
            a1fminus=a1fminus+ccminus(k)*ximinus(1,k)
     +           *cexp(lambdaminus(k)*cmplx(z,0.0))
            a1bminus=a1bminus+ccminus(k)*ximinus(2,k)
     +           *cexp(lambdaminus(k)*cmplx(z,0.0))
            a2fminus=a2fminus+ccminus(k)*ximinus(3,k)
     +           *cexp(lambdaminus(k)*cmplx(z,0.0))
 50      continue
      endif

      return
      end
