c
c
c
      complex*16 function mytanh(z)
      complex*16 z

      if (abs(real(real(z))) .lt. 25) then
         mytanh=(exp(z)-exp(-z))/(exp(z)+exp(-z))
      elseif (real(real(z)) .lt. 25) then
         mytanh=-1.0
      else
         mytanh=1.0
      endif
      return
      end
c
c
c

      subroutine mexFunction(nlhs,plhs,nrhs,prhs)
      integer plhs(*),prhs(*)
      integer nlhs,nrhs
c
c
c
      integer mdim,ndim,xpr,xpi
      integer mxr,mxdelta,mxh,mxm,mxmu,mxsigma,mxsigma0,mxmu0,mxd
      integer mxomega
      real*8 temp
c
c Variables in which we do the actual computations.
c
      integer m
      real*8  lambda(1000),R,DELTA,H,MU0,MU(1000),SIGMA(1000)
      real*8  SIGMA0,D(1000)
      real*8  OMEGA
      complex*16 N(1000),U(1000),Y(1000),N0,R0(1000)
      complex*16 imagunit
c
c Functions
c
      complex*16 mytanh
      integer mxGetM,mxGetN,mxGetPr,mxGetPi
      integer  mxCreateFull

c
c Set up a variable containing i for future use.
c
      imagunit=cmplx(0.0,1.0)
c
c Check the arguments.
c
      if (nlhs .ne. 1) then
         call mexErrMsgTxt('One input required.')
      endif
      if (nrhs .ne. 11) then
         call mexErrMsgTxt('One output required.')
      endif

c
c The first input is a vector of lambdas.
c
      mdim=mxGetM(prhs(1))
      ndim=mxGetN(prhs(1))
      
      if ((mdim .lt. 1) .OR. (ndim .NE. 1)) then
         call mexErrMsgTxt('Input must be a column vector')
      endif

      if (mdim .gt. 1000) then
         call mexErrMsgTxt('Column vector must have < 1000 elements')
      endif

      if (mxIsNumeric(prhs(1)) .eq. 0) then
         call mexErrMsgTxt('Input must be a number')
      endif

c
c Copy the first argument into the subroutine.
c
      xpr=mxGetPr(prhs(1))
      call mxCopyPtrToReal8(xpr,lambda,mdim)

c
c Get the remaining inputs.
c
      mxm=mxGetPr(prhs(5))
      call mxCopyPtrToReal8(mxm,temp,1)
      M=temp

      mxr=mxGetPr(prhs(2))
      call mxCopyPtrToReal8(mxr,R,1)

      mxdelta=mxGetPr(prhs(3))
      call mxCopyPtrToReal8(mxdelta,DELTA,1)

      mxh=mxGetPr(prhs(4))
      call mxCopyPtrToReal8(mxdelta,H,1)

      mxmu=mxGetPr(prhs(6))
      call mxCopyPtrToReal8(mxmu,MU,M)

      mxmu0=mxGetPr(prhs(7))
      call mxCopyPtrToReal8(mxmu0,MU0,1)

      mxsigma=mxGetPr(prhs(8))
      call mxCopyPtrToReal8(mxsigma,SIGMA,M)

      mxsigma0=mxGetPr(prhs(9))
      call mxCopyPtrToReal8(mxsigma0,SIGMA0,1)

      mxd=mxGetPr(prhs(10))
      call mxCopyPtrToReal8(mxd,D,M-1)

      mxomega=mxGetPr(prhs(11))
      call mxCopyPtrToReal8(mxomega,OMEGA,1)

c
c Create a space for the returned value.
c
      plhs(1)=mxCreateFull(mdim,1,1)


c
c Now, do the actual calculation, looping through lambdas.
c
      DO 30 I=1,mdim
c
c First, calculate N and u values.
c
      DO 10 k=1,M
         N(k)=sqrt(lambda(I)**2+imagunit*SIGMA(k)*MU(k)*OMEGA)/
     1            (imagunit*MU(K)*OMEGA)
         u(k)=sqrt(lambda(I)**2+imagunit*SIGMA(k)*MU(k)*OMEGA)
 10   CONTINUE

c
c Then, calculate the Y's.
c
      Y(M-1)=N(M-1)*(N(M)+N(M-1)*mytanh(u(M-1)*D(M-1)))/
     1       (N(M-1)+N(M)*mytanh(u(M-1)*D(M-1)))
      DO 20 K=M-2,1,-1
         Y(k)=N(k)*(Y(k+1)+N(k)*mytanh(u(k)*D(k)))/
     1        (N(k)+Y(k+1)*mytanh(u(k)*D(k)))
 20   CONTINUE

c
c Next, N0
c
      N0=sqrt(lambda(I)**2+imagunit*SIGMA0*MU0*OMEGA)/
     1   (imagunit*MU0*OMEGA)
c
c Finally, the result.
c
      R0(I)=(N0-Y(1))/(N0+Y(1))

c
c End of main loop.
c
 30   CONTINUE
c
c Copy the result back into the returned value.
c
      xpr=mxGetPr(plhs(1))
      xpi=mxGetPi(plhs(1))
      call mxCopyComplex16ToPtr(R0,xpr,xpi,mdim)
c
c Get out of here...
c
      return
      end

