C ALGORITHM 775, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 23,NO. 4, December, 1997, P. 453--493.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# Driver
# Info
# Src
# This archive created: Fri Jul 31 08:51:19 1998
export PATH; PATH=/bin:$PATH
if test ! -d 'Driver'
then
mkdir 'Driver'
fi
cd 'Driver'
if test ! -d 'Dp'
then
mkdir 'Dp'
fi
cd 'Dp'
if test -f 'driver.f'
then
echo shar: will not over-write existing file "'driver.f'"
else
cat << \SHAR_EOF > 'driver.f'
PROGRAM MAIN
C ****************************************************************************
C *
C IMPORTANT NOTES *
C *************** *
C *
C 1. IWORK doesn't have to be as big as it is set here unless you are *
C running at seriously tight tolerances. *
C *
C 2. CPU times are measured using the routine DTIME, which is available *
C on SUN and Silicon Graphics machines (although it doesn't always *
C seem to do what the documentation says it does). For other machines, *
C DTIME will have to be replaced. *
C *
C 3. SLEUTH itself does not need the array WORKS to be as big as it is. *
C It needs an array of length 1136. The big size of WORKS is due to *
C the computation of eigenfunctions later on, in the call to SL4EFN. *
C SL4EFN requires 27 DOUBLE PRECISION storage allocations for each *
C meshpoint used, plus a further 56 DOUBLE PRECISION storage allocations. *
C We could have written SL4EFN to require very few storage allocations *
C but then the eigenfunction computation process would either have been *
C more expensive or else less flexible in the way that the user can *
C specify the evaluation points. *
C *
C *
C END IMPORTANT NOTES *
C ****************************************************************************
C .. Parameters ..
C
INTEGER IWORK
PARAMETER (IWORK=40960)
INTEGER IWORKS
PARAMETER (IWORKS=27*IWORK+56)
C ..
C .. Scalars in Common ..
DOUBLE PRECISION ALFA,BETA,M,RPARM
INTEGER IPROB
C ..
C .. Local Scalars ..
DOUBLE PRECISION A,B,ELAM,ERR,TOL,X
REAL ELAPSE
INTEGER IFAIL,IK,IMATCH,K,MULTI,NMESH,NXTRAP
LOGICAL FSTPAS,HOT,SYM
C ..
C .. Local Arrays ..
DOUBLE PRECISION AEND(10),BEND(10),TARRAY(2),WORK(0:IWORK,1:6),
& WORKS(IWORKS),Y1(1:4),Y2(1:4),YSTOR1(4,0:IWORK),
& YSTOR2(4,0:IWORK)
INTEGER KINDEX(10,10),NUMEIG(10)
LOGICAL SYMARR(10)
C ..
C .. External Functions ..
REAL DTIME
EXTERNAL DTIME
C ..
C .. External Subroutines ..
EXTERNAL SL4BCS,SL4COF,SL4EFN,SLEUTH
C ..
C .. Intrinsic Functions ..
INTRINSIC ATAN
C ..
C .. Common blocks ..
COMMON /COMBLK/IPROB
COMMON /PARM/ALFA,RPARM,M,BETA
C ..
C This code runs the test-problems described in the paper
C `The Code SLEUTH for solving Fourth Order Sturm-Liouville
C Problems', by Greenberg and Marletta. All the input is
C hardwired (see below). The output for each test is as
C follows:
C
C ** The number of the problem attempted
C ** The index of the eigenvalue sought
C ** The error flag
C ** The eigenvalue approximation and error assessment
C ** The CPU time required
C ** The number of Richardson extrapolations used
C ** The number of mesh intervals used
C ** The value of the eigenfunction at the centre of the
C interval [A,B] on which the problem is posed.
C
C Arrays of information about the problems (endpoints,
C number of eigenvalues to be computed and their
C indices):
C For problem 1:
AEND(1) = -2.D0*ATAN(1.D0)
BEND(1) = -AEND(1)
SYMARR(1) = .TRUE.
NUMEIG(1) = 3
KINDEX(1,1) = 2
KINDEX(1,2) = 3
KINDEX(1,3) = 4
C For problem 2:
AEND(2) = 0.D0
BEND(2) = 4.D0*ATAN(1.D0)
SYMARR(2) = .FALSE.
NUMEIG(2) = 3
KINDEX(2,1) = 0
KINDEX(2,2) = 9
KINDEX(2,3) = 99
C For problem 3:
AEND(3) = -100.D0
BEND(3) = 100.D0
SYMARR(3) = .TRUE.
NUMEIG(3) = 3
KINDEX(3,1) = 0
KINDEX(3,2) = 9
KINDEX(3,3) = 99
C For problem 4:
AEND(4) = 0.D0
BEND(4) = 100.D0
SYMARR(4) = .FALSE.
NUMEIG(4) = 3
KINDEX(4,1) = 0
KINDEX(4,2) = 9
KINDEX(4,3) = 19
C For problem 5:
AEND(5) = -0.999999D0
BEND(5) = 0.999999D0
SYMARR(5) = .TRUE.
NUMEIG(5) = 5
KINDEX(5,1) = 0
KINDEX(5,2) = 1
KINDEX(5,3) = 2
KINDEX(5,4) = 3
KINDEX(5,5) = 9
C For problem 6:
AEND(6) = 1.0D-7
BEND(6) = 1.D0
SYMARR(6) = .FALSE.
NUMEIG(6) = 5
KINDEX(6,1) = 0
KINDEX(6,2) = 1
KINDEX(6,3) = 2
KINDEX(6,4) = 3
KINDEX(6,5) = 9
C For problem 7:
AEND(7) = -0.999999D0
BEND(7) = 0.999999D0
SYMARR(7) = .TRUE.
NUMEIG(7) = 5
KINDEX(7,1) = 0
KINDEX(7,2) = 1
KINDEX(7,3) = 2
KINDEX(7,4) = 3
KINDEX(7,5) = 9
C For problem 8:
AEND(8) = 1.0D-6
BEND(8) = 30.D0
SYMARR(8) = .FALSE.
NUMEIG(8) = 5
KINDEX(8,1) = 0
KINDEX(8,2) = 1
KINDEX(8,3) = 2
KINDEX(8,4) = 3
KINDEX(8,5) = 9
C For problem 9:
AEND(9) = 1.0D-6
BEND(9) = 0.999999D0
SYMARR(9) = .FALSE.
NUMEIG(9) = 5
KINDEX(9,1) = 0
KINDEX(9,2) = 1
KINDEX(9,3) = 2
KINDEX(9,4) = 3
KINDEX(9,5) = 9
C For problem 10:
AEND(10) = -0.999999D0
BEND(10) = 0.999999D0
SYMARR(10) = .TRUE.
NUMEIG(10) = 5
KINDEX(10,1) = 0
KINDEX(10,2) = 1
KINDEX(10,3) = 2
KINDEX(10,4) = 3
KINDEX(10,5) = 9
C Tolerance:
TOL = 1.0D-6
FSTPAS = .TRUE.
DO 30 IPROB = 1,10
WRITE (6,FMT=*)
WRITE (7,FMT=*)
WRITE (6,FMT=*) 'Problem Number:',IPROB
WRITE (6,FMT=*) '======================================'
WRITE (7,FMT=*) 'Problem Number:',IPROB
WRITE (7,FMT=*) '======================================'
10 IF (IPROB.EQ.1) THEN
C Problem 1 is done with two different values for the parameter
C BETA:
IF (FSTPAS) THEN
BETA = 10.D0
ELSE
BETA = 30.D0
END IF
END IF
IF (IPROB.EQ.7) THEN
ALFA = 1.D0
END IF
IF (IPROB.EQ.8) THEN
RPARM = 1.D0
END IF
IF (IPROB.EQ.9) THEN
ALFA = 1.D0
M = 1.D0
END IF
C For each problem, compute the eigenvalues with the indices
C K specified:
DO 20 IK = 1,NUMEIG(IPROB)
K = KINDEX(IPROB,IK)
WRITE (6,FMT=*) 'Eigenvalue index:',K
WRITE (6,FMT=*) '-------------------------------------'
WRITE (7,FMT=*) 'Eigenvalue index:',K
WRITE (7,FMT=*) '-------------------------------------'
A = AEND(IPROB)
B = BEND(IPROB)
SYM = SYMARR(IPROB)
ELAPSE = DTIME(TARRAY)
C
C Pretty random-looking initial guess for the eigenvalue;
C did not use ELAM = 0 since some of the problems have
C lambda = 0 as an eigenvalue, which would make life too
C easy for the code.
C
ELAM = 0.123456D0
ERR = 1.D-1
IFAIL = 0
CALL SLEUTH(A,B,ELAM,ERR,K,SYM,SL4COF,SL4BCS,TOL,NMESH,
& IMATCH,NXTRAP,WORK,IWORK,WORKS,IFAIL)
ELAPSE = DTIME(TARRAY) - ELAPSE
IF (IFAIL.EQ.0.D0) THEN
WRITE (6,FMT=*) 'Successful exit from SLEUTH'
WRITE (6,FMT=*) 'Eigenvalue approximation:',ELAM
WRITE (6,FMT=*) 'Error estimate:',ERR
WRITE (6,FMT=*) 'Number of extrapolations:',NXTRAP
WRITE (6,FMT=*) 'CPU Time:',ELAPSE
WRITE (6,FMT=*) 'Number of meshpoints:',NMESH
WRITE (6,FMT=*) 'Matchpoint:',IMATCH
WRITE (7,FMT=*) 'Successful exit from SLEUTH'
WRITE (7,FMT=*) 'Eigenvalue approximation:',ELAM
WRITE (7,FMT=*) 'Error estimate:',ERR
WRITE (7,FMT=*) 'Number of extrapolations:',NXTRAP
WRITE (7,FMT=*) 'CPU Time:',ELAPSE
WRITE (7,FMT=*) 'Number of meshpoints:',NMESH
WRITE (7,FMT=*) 'Matchpoint:',IMATCH
ELSE
WRITE (6,FMT=*) 'Failure exit from SLEUTH, IFAIL=',
& IFAIL
WRITE (6,FMT=*) 'Eigenvalue approximation:',ELAM
WRITE (6,FMT=*) 'Error estimate:',ERR
WRITE (6,FMT=*) 'Number of extrapolations:',NXTRAP
WRITE (6,FMT=*) 'CPU Time:',ELAPSE
WRITE (7,FMT=*) 'Failure exit from SLEUTH, IFAIL=',
& IFAIL
WRITE (7,FMT=*) 'Eigenvalue approximation:',ELAM
WRITE (7,FMT=*) 'Error estimate:',ERR
WRITE (7,FMT=*) 'Number of extrapolations:',NXTRAP
WRITE (7,FMT=*) 'CPU Time:',ELAPSE
END IF
IF (IFAIL.EQ.0) THEN
C Compute the eigenfunction at the midpoint:
HOT = .FALSE.
X = (WORK(0,1)+WORK(NMESH,1))*0.5D0
CALL SL4EFN(X,Y1,Y2,HOT,WORK,IWORK,NMESH,IMATCH,ELAM,
& MULTI,WORKS,YSTOR1,YSTOR2,SL4BCS,IFAIL)
IF (IFAIL.NE.0) THEN
WRITE (6,FMT=*)
& 'Failure exit from SL4EFN, IFAIL =',IFAIL
WRITE (7,FMT=*)
& 'Failure exit from SL4EFN, IFAIL =',IFAIL
END IF
WRITE (6,FMT=*) 'Multiplicity was',MULTI
WRITE (7,FMT=*) 'Multiplicity was',MULTI
WRITE (6,FMT=*) X,Y1(1)
WRITE (7,FMT=*) X,Y1(1)
END IF
20 CONTINUE
IF (IPROB.EQ.1 .AND. FSTPAS) THEN
FSTPAS = .FALSE.
GO TO 10
END IF
30 CONTINUE
STOP
END
C -------------------------------------------------------------------------
SUBROUTINE SL4COF(X,P,S,Q,W)
C .. Scalars in Common ..
DOUBLE PRECISION ALFA,BETA,M,RPARM
INTEGER IPROB
C ..
C .. Common blocks ..
COMMON /COMBLK/IPROB
COMMON /PARM/ALFA,RPARM,M,BETA
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION P,Q,S,W,X
C ..
C .. Intrinsic Functions ..
INTRINSIC COS,EXP,SIN
C ..
GO TO (10,20,30,40,50,60,70,80,90,
& 100,110,120,130) IPROB
C Coffey Evans squared
10 W = 1.D0
P = 1.D0
S = 2.D0* ((BETA*SIN(2.D0*X))**2) - 4.D0*BETA*COS(2.D0*X)
Q = ((BETA**2)* (SIN(2.D0*X)**2)-2.D0*BETA*COS(2.D0*X))**2 -
& 8.D0* (BETA**2)* ((COS(2.D0*X)**2)- (SIN(2.D0*X)**2)) -
& 8.D0*BETA*COS(2.D0*X)
RETURN
C Made-up-problem (oscillatory q(x)) squared
20 P = 1.D0
W = 1.D0
S = 2.D0*COS(X) + 4.D0*COS(2.D0*X) + 6.D0*COS(3.D0*X)
Q = (COS(X)+2.D0*COS(2.D0*X)+3.D0*COS(3.D0*X))**2 + COS(X) +
& 8.D0*COS(2.D0*X) + 27.D0*COS(3.D0*X)
RETURN
C Modified harmonic oscillator squared
30 P = 1.D0
W = 1.D0
S = 2.D0* ((X**2)+ (X**4))
Q = ((X**2)+ (X**4))**2 - 2.D0 - 12.D0* (X**2)
RETURN
C Harmonic oscillator squared
40 P = 1.D0
W = 1.D0
S = 2.D0* (X**2)
Q = X**4 - 2
RETURN
C Legendre equation squared
50 P = (1-X**2.D0)**2.D0
W = 1.D0
Q = 0.D0
S = 2.D0* (1-X**2.D0)
RETURN
C Bessel equation squared
60 P = X
S = 9.D0/X
Q = 0.D0
W = X
RETURN
C Littlejohn Problem 1
70 W = 1.D0
P = 1.D0 - X*X
P = P*P
S = 4.D0*ALFA* (1.D0-X*X) + 8.D0
Q = 0.D0
RETURN
C Littlejohn Problem 2
80 W = EXP(-X)
P = X*X*W
S = ((2.D0*RPARM+2.D0)*X+2.D0)*W
Q = 0.D0
RETURN
C Littlejohn Problem 3
90 P = X*X* ((1.D0-X)** (2.D0+ALFA))
S = 2.D0* ((1.D0-X)** (1.D0+ALFA))* ((ALFA+M+1.D0)*X+1.D0)
Q = 0.D0
W = (1.D0-X)**ALFA
RETURN
C Littlejohn Problem 1 with ALFA = 0
100 W = 1.D0
P = 1.D0 - X*X
P = P*P
S = 8.D0
Q = 0.D0
RETURN
C Additional Problem 1
110 P = 1.D0
S = X + X
Q = X*X
W = 1.D0
RETURN
C Additional Problem 2
120 P = 1.D0
S = X + X
Q = X*X
W = 1.D0
RETURN
C Additional Problem 3
130 P = 1.D0
S = 0.D0
Q = X*X
W = 1.D0
RETURN
END
SUBROUTINE SL4BCS(IEND,ISING,X,U,V,ELAM)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,X
INTEGER IEND
LOGICAL ISING
C ..
C .. Array Arguments ..
DOUBLE PRECISION U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION ALPHA
C ..
C .. Scalars in Common ..
DOUBLE PRECISION ALFA,BETA,M,RPARM
INTEGER IPROB
C ..
C .. Common blocks ..
COMMON /COMBLK/IPROB
COMMON /PARM/ALFA,RPARM,M,BETA
C ..
C .. Intrinsic Functions ..
INTRINSIC EXP,SIN
C ..
GO TO (10,20,30,40,50,60,70,80,90,
& 100,110,120,130) IPROB
10 U(1,1) = 0.D0
U(1,2) = 0.D0
U(2,1) = 1.D0
U(2,2) = 0.D0
V(1,1) = 1.D0
V(1,2) = 1.D0
V(2,1) = 0.D0
V(2,2) = 0.D0
RETURN
20 IF (IEND.EQ.0) THEN
U(1,1) = 0.D0
U(1,2) = 0.D0
U(2,1) = 1.D0
U(2,2) = 0.D0
V(1,1) = 1.D0
V(1,2) = 1.D0
V(2,1) = 0.D0
V(2,2) = 0.D0
ELSE
ALPHA = SIN(X) + 4.D0*SIN(2.D0*X) + 9.D0*SIN(3.D0*X)
U(1,1) = 0.D0
U(1,2) = 1.D0
U(2,1) = 0.D0
U(2,2) = 0.D0
V(1,1) = 0.D0
V(1,2) = ALPHA
V(2,1) = 1.D0
V(2,2) = 0.D0
END IF
RETURN
30 U(1,1) = 0.D0
U(1,2) = 0.D0
U(2,1) = 1.D0
U(2,2) = 0.D0
V(1,1) = 1.D0
V(1,2) = 1.D0
V(2,1) = 0.D0
V(2,2) = 0.D0
RETURN
C40 U(1,1) = 0.D0
C U(1,2) = 0.D0
C U(2,1) = 1.D0
C U(2,2) = 0.D0
C V(1,1) = 1.D0
C V(1,2) = 1.D0
C V(2,1) = 0.D0
C V(2,2) = 0.D0
C RETURN
40 U(1,1) = 1.D0
U(1,2) = 0.D0
U(2,1) = 0.D0
U(2,2) = 1.D0
V(1,1) = 0.D0
V(1,2) = 0.D0
V(2,1) = 0.D0
V(2,2) = 0.D0
RETURN
50 U(1,1) = 1.D0
U(1,2) = 0.D0
U(2,1) = 0.D0
U(2,2) = 0.D0
V(1,1) = 0.D0
V(1,2) = 0.D0
V(2,1) = 0.D0
V(2,2) = 1.D0
RETURN
60 IF (IEND.EQ.0) THEN
ALPHA = (1.D0/X)* (2.D0+ELAM*X*X)/ (1.D0+ (ELAM/4.D0)*X**2.D0)
C WRITE (10,FMT=*) 'EP = ',X,' AND ALPHA = ',ALPHA
C ALPHA = 2.D0/X
U(1,1) = 1.D0
U(1,2) = 0.D0
U(2,1) = ALPHA
U(2,2) = 0.D0
V(1,1) = 0.D0
V(1,2) = 1.D0
V(2,1) = (1.D0/ALPHA)* ((8.D0/X)+8.D0*ALPHA- (ALPHA**2.D0))
V(2,2) = -1.D0/ALPHA
ELSE
U(1,1) = 0.D0
U(1,2) = 0.D0
U(2,1) = 1.D0
U(2,2) = 0.D0
V(1,1) = 0.D0
V(1,2) = 1.D0
V(2,1) = -1.D0
V(2,2) = 0.D0
END IF
RETURN
C Lambda-dependent BC for Littlejohn Problem 1
70 U(1,1) = 1.D0
U(2,1) = ELAM/ (8.D0*ALFA*X)
U(1,2) = 0.D0
U(2,2) = 0.D0
V(1,1) = ELAM/ (ALFA*X)
V(2,1) = 0.D0
V(1,2) = -ELAM/ (8.D0*ALFA*X)
V(2,2) = 1.D0
RETURN
C Lambda-dependent BC for Littlejohn Problem 2
80 IF (IEND.EQ.0) THEN
C U(1,1) = 1.D0
C U(2,1) = -0.5D0*ELAM/RPARM
C U(1,2) = 0.D0
C U(2,2) = 0.D0
C V(1,1) = -ELAM/RPARM
C V(2,1) = 0.D0
C V(1,2) = 0.5D0*ELAM/RPARM
C V(2,2) = 1.D0
U(1,1) = 1.D0
U(2,1) = 0.D0
U(1,2) = -3.D0 + 3.D0*X*X/10.D0
U(2,2) = 3.D0*X/10.D0
V(1,1) = 0.D0
V(2,1) = 0.D0
V(1,2) = 3.D0*X*X*EXP(-X)
V(2,2) = 3.D0*X*X*EXP(-X)
ELSE
U(1,1) = 1.D0
U(1,2) = 0.D0
U(2,1) = 0.D0
U(2,2) = 1.D0
V(1,1) = 0.D0
V(2,1) = 0.D0
V(1,2) = 0.D0
V(2,2) = 0.D0
END IF
RETURN
C Lambda-dependent BC for Littlejohn Problem 3
90 IF (IEND.EQ.0) THEN
U(1,1) = 1.D0
U(2,1) = -0.5D0*ELAM/M
U(1,2) = 0.D0
U(2,2) = 0.D0
V(1,1) = -ELAM/M
V(2,1) = 0.D0
V(1,2) = 0.5D0*ELAM/M
V(2,2) = 1.D0
ELSE
U(1,1) = 1.D0
U(2,1) = 0.D0
U(1,2) = 0.D0
U(2,2) = 1.D0
V(1,1) = 0.D0
V(2,1) = 0.D0
V(1,2) = 0.D0
V(2,2) = 0.D0
END IF
RETURN
C Lambda-independent BC for Littlejohn Problem 1 with ALFA = 0.
100 U(1,1) = 0.D0
U(2,1) = 1.D0
U(1,2) = 0.D0
U(2,2) = 0.D0
V(1,1) = 0.D0
V(2,1) = 0.D0
V(1,2) = 1.D0
V(2,2) = 0.D0
RETURN
C Additional problem 1
110 U(1,1) = 1.D0
U(1,2) = 0.D0
U(2,1) = 0.D0
U(2,2) = 0.D0
V(1,1) = -1.D0
V(2,1) = 0.D0
V(1,2) = 0.D0
V(2,2) = 1.D0
RETURN
C Additional problem 2
120 U(1,1) = 0.D0
U(1,2) = 0.D0
U(2,1) = 0.D0
U(2,2) = 0.D0
V(1,1) = 1.D0
V(2,1) = 0.D0
V(1,2) = 0.D0
V(2,2) = 1.D0
RETURN
C Additional problem 3
130 U(1,1) = 0.D0
U(1,2) = 0.D0
U(2,1) = 0.D0
U(2,2) = 0.D0
V(1,1) = 1.D0
V(2,1) = 0.D0
V(1,2) = 0.D0
V(2,2) = 1.D0
RETURN
END
SHAR_EOF
chmod +x 'driver.f'
fi # end of overwriting check
if test -f 'RES'
then
echo shar: will not over-write existing file "'RES'"
else
cat << \SHAR_EOF > 'RES'
Problem Number: 1
======================================
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 56 4879.7642088549
Successful exit from SLEUTH
Eigenvalue approximation: 4871.3813012287
Error estimate: 4.0841833106242D-04
Number of extrapolations: 4
CPU Time: 1.86658
Number of meshpoints: 448
Matchpoint: 224
Multiplicity was 1
-1.1102230246252D-16 0.68041212746720
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 60 4977.9433464156
Successful exit from SLEUTH
Eigenvalue approximation: 4976.9512353526
Error estimate: 1.5226044967979D-03
Number of extrapolations: 3
CPU Time: 1.43839
Number of meshpoints: 240
Matchpoint: 120
Multiplicity was 1
-1.1102230246252D-16 -4.1467110718111D-09
Eigenvalue index: 4
-------------------------------------
NMESH,ELAM1: 58 5106.8040644204
Successful exit from SLEUTH
Eigenvalue approximation: 5098.7101275991
Error estimate: 4.4682584314311D-03
Number of extrapolations: 3
CPU Time: 0.923593
Number of meshpoints: 232
Matchpoint: 116
Multiplicity was 1
-1.1102230246252D-16 0.72410999798496
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 120 53697.038234892
Successful exit from SLEUTH
Eigenvalue approximation: 53668.639538172
Error estimate: 2.6718120245884D-03
Number of extrapolations: 3
CPU Time: 1.73744
Number of meshpoints: 480
Matchpoint: 240
Multiplicity was 1
-1.1102230246252D-16 1.4346868883653
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 142 53742.780754913
NMESH,ELAM1: 284 53687.378557614
Successful exit from SLEUTH
Eigenvalue approximation: 53668.639438278
Error estimate: 1.0577574512960D-03
Number of extrapolations: 3
CPU Time: 15.5696
Number of meshpoints: 1136
Matchpoint: 568
Multiplicity was 1
-1.1102230246252D-16 0.83696351470774
Eigenvalue index: 4
-------------------------------------
NMESH,ELAM1: 156 53719.470031326
NMESH,ELAM1: 312 53681.343025394
Successful exit from SLEUTH
Eigenvalue approximation: 53668.639473767
Error estimate: 2.0622485802354D-05
Number of extrapolations: 3
CPU Time: 4.43958
Number of meshpoints: 1248
Matchpoint: 624
Multiplicity was 1
-1.1102230246252D-16 1.4345067318102
Problem Number: 2
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 40 0.54642905753509
Successful exit from SLEUTH
Eigenvalue approximation: 0.53638795825054
Error estimate: 7.9018456703809D-07
Number of extrapolations: 4
CPU Time: 1.21062
Number of meshpoints: 320
Matchpoint: 156
Multiplicity was 1
1.5707963267949 -0.76997716667125
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 41 8148.9164904450
Successful exit from SLEUTH
Eigenvalue approximation: 8148.6631335641
Error estimate: 2.0179505620642D-03
Number of extrapolations: 3
CPU Time: 0.667296
Number of meshpoints: 164
Matchpoint: 39
Multiplicity was 1
1.5707963267949 0.56553305631384
Eigenvalue index: 99
-------------------------------------
NMESH,ELAM1: 290 98014968.524860
Successful exit from SLEUTH
Eigenvalue approximation: 98014951.931167
Error estimate: 4.1484231650829
Number of extrapolations: 2
CPU Time: 2.16724
Number of meshpoints: 580
Matchpoint: 91
Multiplicity was 1
1.5707963267949 -0.59495785749792
Problem Number: 3
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 40 2.2201031422153
Successful exit from SLEUTH
Eigenvalue approximation: 1.9386428402571
Error estimate: 5.1268903306460D-07
Number of extrapolations: 5
CPU Time: 1.63204
Number of meshpoints: 640
Matchpoint: 320
Multiplicity was 1
0. 0.83473989504501
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 92 2209.6799075483
Successful exit from SLEUTH
Eigenvalue approximation: 2205.7121356824
Error estimate: 2.8526525053773D-04
Number of extrapolations: 4
CPU Time: 2.42855
Number of meshpoints: 736
Matchpoint: 368
Multiplicity was 1
0. -3.9395984851185D-11
Eigenvalue index: 99
-------------------------------------
NMESH,ELAM1: 6990 1044325.8183711
Successful exit from SLEUTH
Eigenvalue approximation: 1044329.6743999
Error estimate: 0.96400719877177
Number of extrapolations: 2
CPU Time: 73.1019
Number of meshpoints: 13980
Matchpoint: 6990
Multiplicity was 1
0. -6.0366317184339D-13
Problem Number: 4
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 40 4.4825040429826D-02
Successful exit from SLEUTH
Eigenvalue approximation: 1.8028428261970D-06
Error estimate: 3.6488338707912D-07
Number of extrapolations: 4
CPU Time: -3.00950
Number of meshpoints: 320
Matchpoint: 1
Multiplicity was 1
49.999999999999 0.
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 65 1297.3256366106
Successful exit from SLEUTH
Eigenvalue approximation: 1296.0000185709
Error estimate: 9.8926073587791D-05
Number of extrapolations: 4
CPU Time: 4.03933
Number of meshpoints: 520
Matchpoint: 1
Multiplicity was 1
50.000000000000 0.
Eigenvalue index: 19
-------------------------------------
NMESH,ELAM1: 109 5776.3735407344
Successful exit from SLEUTH
Eigenvalue approximation: 5775.9998223221
Error estimate: 1.8576065524810D-04
Number of extrapolations: 4
CPU Time: 5.85964
Number of meshpoints: 872
Matchpoint: 1
Multiplicity was 1
50.000000000000 0.
Problem Number: 5
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 40 5.3819872830008D-13
Successful exit from SLEUTH
Eigenvalue approximation: -1.7862867960038D-13
Error estimate: 1.7920685197512D-13
Number of extrapolations: 2
CPU Time: 2.84470E-02
Number of meshpoints: 80
Matchpoint: 40
Multiplicity was 1
1.1102230246252D-16 0.70710713474020
Eigenvalue index: 1
-------------------------------------
NMESH,ELAM1: 144 4.0015257818330
Successful exit from SLEUTH
Eigenvalue approximation: 4.0000240796134
Error estimate: 5.0753040685692D-07
Number of extrapolations: 3
CPU Time: 4.35831
Number of meshpoints: 576
Matchpoint: 288
Multiplicity was 1
1.1102230246252D-16 3.6839080403047D-10
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 184 35.963849457226
Successful exit from SLEUTH
Eigenvalue approximation: 36.000360361573
Error estimate: 5.6623564342810D-06
Number of extrapolations: 3
CPU Time: 1.61496
Number of meshpoints: 736
Matchpoint: 368
Multiplicity was 1
1.1102230246252D-16 0.79193216681046
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 194 143.08557504765
Successful exit from SLEUTH
Eigenvalue approximation: 144.00201571767
Error estimate: 3.4194420737776D-05
Number of extrapolations: 4
CPU Time: 5.19443
Number of meshpoints: 1552
Matchpoint: 776
Multiplicity was 1
1.1102230246252D-16 1.5555413524675D-15
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 344 8093.5449934735
Successful exit from SLEUTH
Eigenvalue approximation: 8100.3076456834
Error estimate: 1.3419749167345D-03
Number of extrapolations: 3
CPU Time: 4.86541
Number of meshpoints: 1376
Matchpoint: 688
Multiplicity was 1
1.1102230246252D-16 -5.3758807620041D-11
Problem Number: 6
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 55 695.21205644539
Successful exit from SLEUTH
Eigenvalue approximation: 695.62040633147
Error estimate: 9.0420013687738D-05
Number of extrapolations: 4
CPU Time: 3.90823
Number of meshpoints: 440
Matchpoint: 439
Multiplicity was 1
0.50000005000000 1.8949018403124
Eigenvalue index: 1
-------------------------------------
NMESH,ELAM1: 54 5016.5906098657
Successful exit from SLEUTH
Eigenvalue approximation: 5019.7261245463
Error estimate: 1.9908896403649D-03
Number of extrapolations: 4
CPU Time: 1.84142
Number of meshpoints: 432
Matchpoint: 431
Multiplicity was 1
0.50000005000000 -1.6059339454802
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 64 18220.604554077
Successful exit from SLEUTH
Eigenvalue approximation: 18230.611376208
Error estimate: 7.3536409106358D-03
Number of extrapolations: 3
CPU Time: 1.02072
Number of meshpoints: 256
Matchpoint: 255
Multiplicity was 1
0.50000005000000 -1.2268382418799
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 80 47912.353555606
Successful exit from SLEUTH
Eigenvalue approximation: 47926.085571259
Error estimate: 9.0248889498374D-03
Number of extrapolations: 3
CPU Time: 1.78353
Number of meshpoints: 320
Matchpoint: 319
Multiplicity was 1
0.50000005000000 1.7094864619280
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 229 1292302.9036365
Successful exit from SLEUTH
Eigenvalue approximation: 1292322.6251748
Error estimate: 4.5412955650439D-03
Number of extrapolations: 3
CPU Time: 3.70746
Number of meshpoints: 916
Matchpoint: 915
Multiplicity was 1
0.50000005000000 -1.7860472098655
Problem Number: 7
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 40 -5.3928733803079D-13
Successful exit from SLEUTH
Eigenvalue approximation: 8.3876841036320D-13
Error estimate: 3.4451393709850D-13
Number of extrapolations: 2
CPU Time: 5.50370E-02
Number of meshpoints: 80
Matchpoint: 40
Multiplicity was 1
1.1102230246252D-16 0.70710713474019
Eigenvalue index: 1
-------------------------------------
NMESH,ELAM1: 98 8.0007587802097
Successful exit from SLEUTH
Eigenvalue approximation: 8.0000060295348
Error estimate: 2.2799505252635D-07
Number of extrapolations: 3
CPU Time: 1.19356
Number of meshpoints: 392
Matchpoint: 196
Multiplicity was 1
1.1102230246252D-16 2.4016240652498D-12
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 136 48.004050650067
Successful exit from SLEUTH
Eigenvalue approximation: 48.000102315643
Error estimate: 4.1240386115267D-06
Number of extrapolations: 3
CPU Time: 1.30555
Number of meshpoints: 544
Matchpoint: 272
Multiplicity was 1
1.1102230246252D-16 1.0542401929870
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 138 167.72599358633
Successful exit from SLEUTH
Eigenvalue approximation: 168.00056027209
Error estimate: 2.5532868594761D-05
Number of extrapolations: 4
CPU Time: 2.28092
Number of meshpoints: 1104
Matchpoint: 552
Multiplicity was 1
1.1102230246252D-16 -5.9567318581469D-12
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 218 8277.9620827909
NMESH,ELAM1: 436 8279.5378489194
Successful exit from SLEUTH
Eigenvalue approximation: 8280.0785838939
Error estimate: 6.2941497890279D-05
Number of extrapolations: 3
CPU Time: 6.83956
Number of meshpoints: 1744
Matchpoint: 872
Multiplicity was 1
1.1102230246252D-16 -6.9157285391312D-14
Problem Number: 8
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 41 -2.7822741905732D-11
Successful exit from SLEUTH
Eigenvalue approximation: -1.1618207715853D-10
Error estimate: 2.2089833813200D-11
Number of extrapolations: 2
CPU Time: 0.970555
Number of meshpoints: 82
Matchpoint: 81
Multiplicity was 1
15.000000500000 1.0008374920648
Eigenvalue index: 1
-------------------------------------
NMESH,ELAM1: 50 5.1301792134227
Successful exit from SLEUTH
Eigenvalue approximation: 5.4160206636243
Error estimate: 8.3834919882366D-07
Number of extrapolations: 5
CPU Time: 6.77207
Number of meshpoints: 800
Matchpoint: 1
Multiplicity was 1
15.000000500000 -22.075130195100
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 65 11.702106417881
NMESH,ELAM1: 130 12.311559852756
Successful exit from SLEUTH
Eigenvalue approximation: 12.475718978329
Error estimate: 1.1408533689779D-05
Number of extrapolations: 4
CPU Time: 9.91049
Number of meshpoints: 1040
Matchpoint: 1
Multiplicity was 1
15.000000500000 125.47395761578
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 66 19.908191314644
NMESH,ELAM1: 132 21.170395585746
Successful exit from SLEUTH
Eigenvalue approximation: 21.498734419417
Error estimate: 1.5945392542941D-06
Number of extrapolations: 5
CPU Time: 25.3963
Number of meshpoints: 2112
Matchpoint: 1
Multiplicity was 1
15.000000500000 -352.29241900810
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 114 134.64958715873
Successful exit from SLEUTH
Eigenvalue approximation: 137.48224274641
Error estimate: 2.5299611237036D-05
Number of extrapolations: 5
CPU Time: 24.4290
Number of meshpoints: 1824
Matchpoint: 1
Multiplicity was 1
15.000000500000 -387.06546358094
Problem Number: 9
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 40 1.4498431861898D-13
Successful exit from SLEUTH
Eigenvalue approximation: 5.9264310128428D-13
Error estimate: 1.1191469566632D-13
Number of extrapolations: 2
CPU Time: -0.497167
Number of meshpoints: 80
Matchpoint: 1
Multiplicity was 1
0.50000000000000 1.4142149765888
Eigenvalue index: 1
-------------------------------------
NMESH,ELAM1: 91 17.915901068056
Successful exit from SLEUTH
Eigenvalue approximation: 18.000034276648
Error estimate: 8.5654949333029D-06
Number of extrapolations: 3
CPU Time: 2.55372
Number of meshpoints: 364
Matchpoint: 1
Multiplicity was 1
0.50000000000000 -1.7271185730746
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 87 86.954920410976
Successful exit from SLEUTH
Eigenvalue approximation: 88.000263555671
Error estimate: 5.2277931599557D-05
Number of extrapolations: 3
CPU Time: 2.48613
Number of meshpoints: 348
Matchpoint: 1
Multiplicity was 1
0.50000000000000 -0.93861568161260
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 96 267.18356344436
Successful exit from SLEUTH
Eigenvalue approximation: 270.00106142964
Error estimate: 1.8628212766695D-04
Number of extrapolations: 4
CPU Time: 4.30589
Number of meshpoints: 768
Matchpoint: 1
Multiplicity was 1
0.50000000000000 1.3851655087114
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 222 10089.9527759388
NMESH,ELAM1: 444 10096.0036208729
Successful exit from SLEUTH
Eigenvalue approximation: 10098.1009837331
Error estimate: 3.0188053133315D-04
Number of extrapolations: 3
CPU Time: 6.98814
Number of meshpoints: 1776
Matchpoint: 1
Multiplicity was 1
0.50000000000000 1.2239777396944
Problem Number: 10
======================================
Eigenvalue index: 0
-------------------------------------
NMESH,ELAM1: 40 23.978590804107
Successful exit from SLEUTH
Eigenvalue approximation: 24.000032402611
Error estimate: 2.0735444390851D-05
Number of extrapolations: 3
CPU Time: -1.30190E-02
Number of meshpoints: 160
Matchpoint: 80
Multiplicity was 1
1.1102230246252D-16 0.96828448758508
Eigenvalue index: 1
-------------------------------------
NMESH,ELAM1: 40 119.25547741589
Successful exit from SLEUTH
Eigenvalue approximation: 120.00039225593
Error estimate: 2.3108918555674D-05
Number of extrapolations: 5
CPU Time: 1.44649
Number of meshpoints: 640
Matchpoint: 320
Multiplicity was 1
1.1102230246252D-16 -7.4693707632871D-12
Eigenvalue index: 2
-------------------------------------
NMESH,ELAM1: 42 356.99593635185
Successful exit from SLEUTH
Eigenvalue approximation: 360.00160980503
Error estimate: 4.1838612401079D-05
Number of extrapolations: 5
CPU Time: 1.50828
Number of meshpoints: 672
Matchpoint: 336
Multiplicity was 1
1.1102230246252D-16 0.84187927575495
Eigenvalue index: 3
-------------------------------------
NMESH,ELAM1: 46 835.46254933378
Successful exit from SLEUTH
Eigenvalue approximation: 840.00461059928
Error estimate: 1.1525974123288D-04
Number of extrapolations: 5
CPU Time: 2.19892
Number of meshpoints: 736
Matchpoint: 368
Multiplicity was 1
1.1102230246252D-16 -5.6752807133559D-15
Eigenvalue index: 9
-------------------------------------
NMESH,ELAM1: 124 17151.496929329
Successful exit from SLEUTH
Eigenvalue approximation: 17160.196917773
Error estimate: 3.4889330015479D-03
Number of extrapolations: 3
CPU Time: 1.84058
Number of meshpoints: 496
Matchpoint: 248
Multiplicity was 1
1.1102230246252D-16 -1.5507586250831D-11
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
if test ! -d 'Info'
then
mkdir 'Info'
fi
cd 'Info'
if test -f 'depend'
then
echo shar: will not over-write existing file "'depend'"
else
cat << \SHAR_EOF > 'depend'
nag/f06yaf
nag/x02ajf
nag/c05azf
nag/f06qhf
nag/f06paf
nag/p01abf
nag/f06qff
nag/x01aaf
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
mkdir 'Src'
fi
cd 'Src'
if test ! -d 'Dp'
then
mkdir 'Dp'
fi
cd 'Dp'
if test -f 'src.f'
then
echo shar: will not over-write existing file "'src.f'"
else
cat << \SHAR_EOF > 'src.f'
C =====================================================================
C ======================= SUBROUTINES FOR EIGENVALUES =================
C =====================================================================
C First revision: 27th April 1995
SUBROUTINE SLEUTH(A,B,ELAM,ERR,K,SYM,SL4COF,SL4BCS,TOL,NMESH,
& IMATCH,NXTRAP,WORK,IWORK,WORKS,IFAIL)
C .. Parameters ..
INTEGER NCOARS,IELAM3
PARAMETER (NCOARS=40,IELAM3=10)
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION A,B,ELAM,ERR,TOL
INTEGER IFAIL,IMATCH,IWORK,K,NMESH,NXTRAP
LOGICAL SYM
C ..
C .. Array Arguments ..
C REMARK: WORKS must be of length at least 27*NCOARS+56
DOUBLE PRECISION WORK(0:IWORK,1:6),WORKS(27*NCOARS+56)
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4BCS,SL4COF
C ..
C .. Local Scalars ..
DOUBLE PRECISION ALOC,BLOC,CORR1,CORR2,EL,ELAM0,ELAM1,ELAM2,EPS,
& EU,MACHPT,MCHEPS,OSC,OSCL,TEORAT,TOL1,TOLI,XSHIFT
INTEGER COUNT,I,IFO,II,IM1ST,IMESH,IREFIN,IY1,IY2,J,JMAX,KNTMAX,
& KR,KXTRAP,MAXTRP,MULTI,NM1ST,NMAX,NREC
LOGICAL RXTRAP
CHARACTER*6 SRNAME
C ..
C .. Local Arrays ..
DOUBLE PRECISION ELAM3(IELAM3,IELAM3)
CHARACTER*80 REC(2)
C ..
C .. External Functions ..
DOUBLE PRECISION X02AJF
INTEGER P01ABF
EXTERNAL X02AJF,P01ABF
C ..
C .. External Subroutines ..
EXTERNAL COARS4,GETVEC,MESH4,SOL4
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,MAX,MIN,SQRT
C ..
IFO = IFAIL
MCHEPS = X02AJF(1.D0)
KXTRAP = IELAM3
MAXTRP = IELAM3
C Error Trapping
IF (TOL.LE.MCHEPS**0.6D0 .OR. A.GE.B .OR. ABS(IFAIL).GT.1 .OR.
& ABS(ERR).LE.0.D0 .OR. K.LT.0 .OR. IWORK.LT.NCOARS) THEN
C We have an input error.
IF (ABS(IFAIL).GT.1) THEN
WRITE (REC,FMT=9000) IFAIL
9000 FORMAT (' ** Illegal initial IFAIL value of',I8,/,
& ' ** IFAIL reset to 1')
IFAIL = 1
IFO = 1
NREC = 2
GO TO 180
END IF
IFAIL = 1
IF (TOL.LE.MCHEPS**0.6D0) THEN
WRITE (REC,FMT=9010) MCHEPS**0.6D0,TOL
9010 FORMAT (' ** On entry TOL was',G18.8,/,
& ' ** TOL must be greater than',G18.8)
NREC = 2
GO TO 180
END IF
IF (A.GE.B) THEN
WRITE (REC,FMT=9020) A,B
9020 FORMAT (' ** A must be less than B; you had',/,' ** A=',
& G18.8,' and B=',G18.8)
NREC = 2
GO TO 180
END IF
IF (ERR.EQ.0.D0) THEN
WRITE (REC,FMT=9030)
9030 FORMAT (' ** ERR must be non-zero on entry')
NREC = 1
GO TO 180
END IF
IF (K.LT.0) THEN
WRITE (REC,FMT=9040) K
9040 FORMAT (' ** K must be non-negative; you had',/,' ** K=',
& I8)
NREC = 2
GO TO 180
END IF
IF (IWORK.LT.20) THEN
WRITE (REC,FMT=9050) IWORK
9050 FORMAT (' ** IWORK must be at least 20;',/,
& ' ** you had IWORK=',I8)
NREC = 2
GO TO 180
END IF
END IF
MCHEPS = X02AJF(0.D0)
C First we get a very crude initial approximation to the eigenvalue
IF (SYM) THEN
XSHIFT = 0.5D0* (A+B)
ELSE
XSHIFT = 0.D0
END IF
NMESH = NCOARS
ALOC = (A-XSHIFT)/ (1.D0+ABS(A-XSHIFT))
BLOC = (B-XSHIFT)/ (1.D0+ABS(B-XSHIFT))
CALL COARS4(NMESH,ALOC,BLOC,WORK(0,2),3)
C Now unmap the mesh back into the original interval
DO 10 I = 0,NMESH
WORK(I,2) = WORK(I,2)/ (1.D0-ABS(WORK(I,2))) + XSHIFT
10 CONTINUE
DO 20 I = 1,NMESH
CALL SL4COF((WORK(I-1,2)+WORK(I,2))/2.D0,WORK(I,3),WORK(I,4),
& WORK(I,5),WORK(I,6))
20 CONTINUE
IMATCH = NMESH/2
ELAM1 = ELAM
EPS = ABS(ERR)
TOL1 = SQRT(MCHEPS)
KNTMAX = 100
IFAIL = 0
CALL SOL4(WORK(0,2),NMESH,IMATCH,WORK(1,3),WORK(1,4),WORK(1,5),
& WORK(1,6),ELAM1,EPS,EL,EU,K,TOL1,KNTMAX,SL4BCS,IFAIL)
IF (IFAIL.NE.0) THEN
IF (IFAIL.EQ.2) THEN
WRITE (REC,FMT=9060)
9060 FORMAT (' ** On first call to SOL4, cannot do zero-count',
& /,' ** correction at centre: unsuitable U and V.')
ELSE IF (IFAIL.EQ.3) THEN
WRITE (REC,FMT=9070)
9070 FORMAT (' ** On first call to SOL4, cannot do zero-count',
& /,
& ' ** correction at centre: overflow computing detU, detV.'
& )
ELSE IF (IFAIL.EQ.4) THEN
WRITE (REC,FMT=9080) IFAIL
9080 FORMAT (
& ' ** On first call to SOL4, an error occurred while'
& ,/,' ** processing a Theta matrix; IFAIL=',I8)
ELSE IF (IFAIL.EQ.6) THEN
WRITE (REC,FMT=9090)
9090 FORMAT (' ** On first call to SOL4, cannot do zero-count',
& /,
& ' ** correction on a mesh interval: unsuitable U and V.'
& )
ELSE IF (IFAIL.EQ.7) THEN
WRITE (REC,FMT=9100)
9100 FORMAT (' ** On first call to SOL4, cannot do zero-count',
& /,' ** correction on a mesh interval: overflow',
& ' computing detU, detV.')
ELSE IF (IFAIL.EQ.9) THEN
WRITE (REC,FMT=9110)
9110 FORMAT (' ** On first call to SOL4, no eigenvalue can be',
& /,' ** found; ERR is too small.')
ELSE IF (IFAIL.EQ.10) THEN
WRITE (REC,FMT=9120)
9120 FORMAT (
& ' ** On first call to SOL4, an eigenvalue has been'
& ,/,
& ' ** bracketed but cannot be accurately located.')
END IF
NREC = 2
ELAM = ELAM1
WORK(0,3) = ELAM1
GO TO 180
END IF
ELAM0 = ELAM1
IY1 = 19* (NMESH+2) + 11
IY2 = IY1 + 4* (NMESH+1)
C WORKS must be of length at least 27*NMESH+65
CALL GETVEC(WORK(0,2),NMESH,IMATCH,WORK(1,3),WORK(1,4),WORK(1,5),
& WORK(1,6),ELAM1,MULTI,WORKS,WORKS(IY1),WORKS(IY2),
& SL4BCS,IFAIL)
C Update the matchpoint
IF (.NOT.SYM) THEN
I = IMATCH
C OSC = (WORK(I,4) + SQRT(MAX(0.D0,WORK(I,4)**2+4.D0*WORK(I,
C & 3)* (ELAM1*WORK(I,6)-WORK(I,5)))))/WORK(I,3)
C DO 30 I = 1,NMESH - 1
C OSCL = (WORK(I,4) + SQRT(MAX(0.D0,WORK(I,4)**2+4.D0*WORK(I,
C & 3)* (ELAM1*WORK(I,6)-WORK(I,5)))))/WORK(I,3)
OSC = 4.D0*WORK(I,3)* (ELAM1*WORK(I,6)-WORK(I,5))/WORK(I,3)
DO 30 I = 1,NMESH - 1
OSCL = 4.D0*WORK(I,3)* (ELAM1*WORK(I,6)-WORK(I,5))/
& WORK(I,3)
IF (OSCL.GT.OSC) THEN
OSC = OSCL
IMATCH = I
END IF
30 CONTINUE
END IF
MACHPT = WORK(IMATCH,2)
C Now use this eigenvalue approximation to assist in meshing.
C Store the matchpoint and the number of mesh-intervals of the coarse
C initial mesh:
IM1ST = IMATCH
NM1ST = NMESH
C Set meshing tolerance for automatic meshing:
C TOLI = 1.D-2
C TOLI = 0.5D0*SQRT(TOL)
C TOLI = TOL
TOLI = 1.D0
COUNT = 0
NMAX = IWORK
40 IFAIL = 0
CALL MESH4(MACHPT,A,ELAM1,SL4COF,WORK(0,1),NMAX,0,IMESH,TOLI,
& WORKS(IY1),WORKS(IY2),MULTI,IM1ST,NM1ST,WORK(0,2),
& IWORK,IFAIL)
IF (IFAIL.NE.0) THEN
IF (COUNT.LT.3) THEN
TOLI = TOLI*5.D0
COUNT = COUNT + 1
GO TO 40
END IF
GO TO 60
END IF
IF (SYM) THEN
C Generate the second half of the mesh by reflection
NMESH = 2*IMESH
DO 50 I = 1,IMESH
WORK(IMESH+I,1) = 2.D0*MACHPT - WORK(IMESH-I,1)
50 CONTINUE
GO TO 80
END IF
CALL MESH4(MACHPT,B,ELAM1,SL4COF,WORK(0,1),NMAX,IMESH,NMESH,TOLI,
& WORKS(IY1),WORKS(IY2),MULTI,IM1ST,NM1ST,WORK(0,2),
& IWORK,IFAIL)
IF (IFAIL.NE.0) THEN
IF (COUNT.LT.3) THEN
TOLI = TOLI*5.D0
COUNT = COUNT + 1
GO TO 40
END IF
END IF
60 IF (IFAIL.EQ.0) THEN
NMESH = NMESH + IMESH
IMATCH = IMESH
ELSE
C The meshing has failed. Adopt contingency measures:
NMESH = NCOARS
IMATCH = NMESH/2
IF (SYM) THEN
XSHIFT = 0.5D0* (A+B)
ELSE
XSHIFT = 0.D0
END IF
ALOC = (A-XSHIFT)/ (1.D0+ABS(A-XSHIFT))
BLOC = (B-XSHIFT)/ (1.D0+ABS(B-XSHIFT))
CALL COARS4(NMESH,ALOC,BLOC,WORK(0,1),3)
C Now unmap the mesh back into the original interval
DO 70 I = 0,NMESH
WORK(I,1) = WORK(I,1)/ (1.D0-ABS(WORK(I,1))) + XSHIFT
70 CONTINUE
END IF
C Evaluate the coefficients
80 DO 90 I = 1,NMESH
CALL SL4COF((WORK(I-1,1)+WORK(I,1))/2.D0,WORK(I,2),WORK(I,3),
& WORK(I,4),WORK(I,5))
90 CONTINUE
IREFIN = 0
RXTRAP = .FALSE.
EPS = ABS(ERR)
C Update the matchpoint again.
100 IF (SYM) THEN
IMATCH = NMESH/2
ELSE
I = IMATCH
C OSC = (WORK(I,3) + SQRT(MAX(0.D0,WORK(I,3)**2+4.D0*WORK(I,
C & 2)* (ELAM1*WORK(I,5)-WORK(I,4)))))/WORK(I,2)
OSC = 4.D0*WORK(I,2)* (ELAM1*WORK(I,5)-WORK(I,4))/WORK(I,2)
DO 110 I = 1,NMESH - 1
C OSCL = (WORK(I,3) + SQRT(MAX(0.D0,WORK(I,3)**2+4.D0*WORK(I,
C & 2)* (ELAM1*WORK(I,5)-WORK(I,4)))))/WORK(I,2)
OSCL = 4.D0*WORK(I,2)* (ELAM1*WORK(I,5)-WORK(I,4))/
& WORK(I,2)
IF (OSCL.GT.OSC) THEN
OSC = OSCL
IMATCH = I
END IF
110 CONTINUE
END IF
C We are now ready to commence computation of decent eigenvalue
C approximations.
IFAIL = 0
TOL1 = MCHEPS**0.75D0
CALL SOL4(WORK(0,1),NMESH,IMATCH,WORK(1,2),WORK(1,3),WORK(1,4),
& WORK(1,5),ELAM1,EPS,EL,EU,K,TOL1,KNTMAX,SL4BCS,IFAIL)
WORK(0,3) = ELAM1
WRITE (6,FMT=*) 'NMESH,ELAM1:',NMESH,ELAM1
IF (IFAIL.NE.0) THEN
IF (IFAIL.EQ.2) THEN
WRITE (REC,FMT=9130)
9130 FORMAT (' ** On 2nd call to SOL4, cannot do zero-count',/,
& ' ** correction at centre: unsuitable U and V.')
ELSE IF (IFAIL.EQ.3) THEN
WRITE (REC,FMT=9140)
9140 FORMAT (' ** On 2nd call to SOL4, cannot do zero-count',/,
& ' ** correction at centre: overflow computing detU, detV.'
& )
ELSE IF (IFAIL.EQ.4) THEN
WRITE (REC,FMT=9150) IFAIL
9150 FORMAT (' ** On 2nd call to SOL4, an error occurred while'
& ,/,' ** processing a Theta matrix; IFAIL=',I8)
ELSE IF (IFAIL.EQ.6) THEN
WRITE (REC,FMT=9160)
9160 FORMAT (' ** On 2nd call to SOL4, cannot do zero-count',/,
& ' ** correction on a mesh interval: unsuitable U and V.'
& )
ELSE IF (IFAIL.EQ.7) THEN
WRITE (REC,FMT=9170)
9170 FORMAT (' ** On 2nd call to SOL4, cannot do zero-count',/,
& ' ** correction on a mesh interval: overflow',
& ' computing detU, detV.')
ELSE IF (IFAIL.EQ.9) THEN
WRITE (REC,FMT=9180)
9180 FORMAT (' ** On 2nd call to SOL4, no eigenvalue can be',/,
& ' ** found; ERR is too small.')
ELSE IF (IFAIL.EQ.10) THEN
WRITE (REC,FMT=9190)
9190 FORMAT (' ** On 2nd call to SOL4, an eigenvalue has been',
& /,' ** bracketed but cannot be accurately located.'
& )
END IF
NREC = 2
ELAM = ELAM1
GO TO 180
END IF
C Update the matchpoint:
IF (.NOT.SYM) THEN
I = IMATCH
C OSC = (WORK(I,3) + SQRT(MAX(0.D0,WORK(I,3)**2+4.D0*WORK(I,
C & 2)* (ELAM1*WORK(I,5)-WORK(I,4)))))/WORK(I,2)
OSC = 4.D0*WORK(I,2)* (ELAM1*WORK(I,5)-WORK(I,4))/WORK(I,2)
DO 120 I = 1,NMESH - 1
C OSCL = (WORK(I,3) + SQRT(MAX(0.D0,WORK(I,3)**2+4.D0*WORK(I,
C & 2)* (ELAM1*WORK(I,5)-WORK(I,4)))))/WORK(I,2)
OSCL = 4.D0*WORK(I,2)* (ELAM1*WORK(I,5)-WORK(I,4))/
& WORK(I,2)
IF (OSCL.GT.OSC) THEN
OSC = OSCL
IMATCH = I
END IF
120 CONTINUE
END IF
EPS = ABS(ELAM1-ELAM0)
ELAM0 = ELAM1
IF (EPS.LT.0.1D0*MAX(1.D0,ABS(ELAM1))) THEN
C The error seems to be small enough to justify commencing
C Richardson extrapolation
ELAM3(1,1) = ELAM1
KR = 2
RXTRAP = .TRUE.
C WRITE (6,FMT=*) 'Started Richardson Extrapolation'
END IF
EPS = MAX(0.1D0,EPS)
130 DO 140 II = NMESH,1,-1
WORK(2*II,1) = WORK(II,1)
WORK(2*II-1,1) = 0.5D0* (WORK(II,1)+WORK(II-1,1))
140 CONTINUE
NMESH = 2*NMESH
DO 150 II = 1,NMESH
CALL SL4COF((WORK(II-1,1)+WORK(II,1))/2.D0,WORK(II,2),
& WORK(II,3),WORK(II,4),WORK(II,5))
150 CONTINUE
C Now update the matchpoint properly:
IF (.NOT.SYM) THEN
I = IMATCH*2
C OSC = (WORK(I,3) + SQRT(MAX(0.D0,WORK(I,3)**2+4.D0*WORK(I,
C & 2)* (ELAM1*WORK(I,5)-WORK(I,4)))))/WORK(I,2)
OSC = 4.D0*WORK(I,2)* (ELAM1*WORK(I,5)-WORK(I,4))/WORK(I,2)
DO 160 I = 1,NMESH - 1
C OSCL = (WORK(I,3) + SQRT(MAX(0.D0,WORK(I,3)**2+4.D0*WORK(I,
C & 2)* (ELAM1*WORK(I,5)-WORK(I,4)))))/WORK(I,2)
OSCL = 4.D0*WORK(I,2)* (ELAM1*WORK(I,5)-WORK(I,4))/
& WORK(I,2)
IF (OSCL.GT.OSC) THEN
OSC = OSCL
IMATCH = I
END IF
160 CONTINUE
ELSE
IMATCH = IMATCH*2
END IF
IF (.NOT.RXTRAP) THEN
IREFIN = IREFIN + 1
IF (IREFIN.LT.10 .AND. 2*NMESH.LT.IWORK) GO TO 100
IFAIL = 8
ELAM = ELAM1
WRITE (REC,FMT=9200)
9200 FORMAT (' ** The eigenvalue approximations are varying',/,
& ' ** wildly as the mesh is refined')
NREC = 2
GO TO 180
END IF
TOL1 = MCHEPS**0.75D0
JMAX = MIN(KXTRAP,KR-1)
ELAM2 = ELAM3(KR-1,JMAX)
EPS = ABS(ERR)
IF (KR.GT.2) EPS = ABS(ELAM3(KR-1,1)-ELAM3(KR-2,1))
IF (KR.GT.2) TOL1 = MAX(TOL1,0.01D0*
& (EPS/ (1.D0+ABS(ELAM3(KR-2,1))))**2)
IFAIL = 0
CALL SOL4(WORK(0,1),NMESH,IMATCH,WORK(1,2),WORK(1,3),WORK(1,4),
& WORK(1,5),ELAM2,EPS,EL,EU,K,TOL1,KNTMAX,SL4BCS,IFAIL)
WORK(0,3) = ELAM2
IF (IFAIL.NE.0) THEN
IF (IFAIL.EQ.2) THEN
WRITE (REC,FMT=9210)
9210 FORMAT (' ** On later call to SOL4, cannot do zero-count',
& /,' ** correction at centre: unsuitable U and V.')
ELSE IF (IFAIL.EQ.3) THEN
WRITE (REC,FMT=9220)
9220 FORMAT (' ** On later call to SOL4, cannot do zero-count',
& /,
& ' ** correction at centre: overflow computing detU, detV.'
& )
ELSE IF (IFAIL.EQ.4) THEN
WRITE (REC,FMT=9230) IFAIL
9230 FORMAT (
& ' ** On later call to SOL4, an error occurred while'
& ,/,' ** processing a Theta matrix; IFAIL=',I8)
ELSE IF (IFAIL.EQ.6) THEN
WRITE (REC,FMT=9240)
9240 FORMAT (' ** On later call to SOL4, cannot do zero-count',
& /,
& ' ** correction on a mesh interval: unsuitable U and V.'
& )
ELSE IF (IFAIL.EQ.7) THEN
WRITE (REC,FMT=9250)
9250 FORMAT (' ** On later call to SOL4, cannot do zero-count',
& /,' ** correction on a mesh interval: overflow',
& 'computing detU, detV.')
ELSE IF (IFAIL.EQ.9) THEN
WRITE (REC,FMT=9260)
9260 FORMAT (' ** On later call to SOL4, no eigenvalue can be',
& /,' ** found; ERR is too small.')
ELSE IF (IFAIL.EQ.10) THEN
WRITE (REC,FMT=9270)
9270 FORMAT (
& ' ** On later call to SOL4, an eigenvalue has been'
& ,/,
& ' ** bracketed but cannot be accurately located.')
END IF
NREC = 2
GO TO 180
END IF
ELAM3(KR,1) = ELAM2
JMAX = MIN(KXTRAP,KR)
C This loop may do one more Richardson extrapolation than is strictly
C justified, but what the hell -- if EPS passes the error test it
C will probably be OK.
DO 170 J = 2,JMAX
EPS = (ELAM3(KR,J-1)-ELAM3(KR-1,J-1))/
& (2.D0** (2* (J-1))-1.D0)
ELAM3(KR,J) = ELAM3(KR,J-1) + EPS
170 CONTINUE
C Check the error:
C ELAM2 = ELAM3(KR,KR)
ELAM2 = ELAM3(KR,JMAX)
EPS = ABS(EPS)
ELAM = ELAM2
ERR = EPS
NXTRAP = KR
IF (EPS.LT.TOL*MAX(1.D0,ABS(ELAM))) THEN
C We have got our approximation
IFAIL = 0
RETURN
END IF
C Now look through the Richardson extrapolation table and see if there
C is evidence of growing Richardson corrections. If the Richardson
C corrections are growing then put an upper limit on the number of
C extrapolations to be attempted.
IF (KR.GT.2) THEN
C First check that the mesh-size is small enough to justify
C extrapolation
CORR2 = ELAM3(KR,1) - ELAM3(KR-1,1)
IF (ABS(CORR2).GT.0.1D0*MAX(1.D0,ABS(ELAM3(KR,1)))) THEN
C WRITE (6,FMT=*)
C & 'The mesh appears too coarse: extrapolation restarted'
C The mesh is too coarse for extrapolation to be justified yet.
IF (KR.GT.MAXTRP) THEN
IFAIL = 8
WRITE (REC,FMT=9280)
9280 FORMAT (
& ' ** More than 10 extrapolations would be required'
& ,/,' ** to obtain the required accuracy')
NREC = 2
GO TO 180
END IF
C ELAM3(1,1) = ELAM3(KR,1)
C MAXTRP = MAXTRP - KR
C KR = 2
C EPS = 2.D0*ABS(CORR2)
RXTRAP = .FALSE.
ELAM1 = ELAM3(1,1)
MAXTRP = MAXTRP - KR
GO TO 130
END IF
END IF
IF (JMAX.GT.2 .AND. KR.GT.2) THEN
CORR2 = (ELAM3(KR,JMAX-2)-ELAM3(KR-1,JMAX-2))/
& MAX(1.D0,ABS(ELAM3(KR,1)))
CORR1 = (ELAM3(KR-1,JMAX-2)-ELAM3(KR-2,JMAX-2))/
& MAX(1.D0,ABS(ELAM3(KR,1)))
TEORAT = 0.5D0** (2* (JMAX-1))
IF (ABS(CORR1).LT.X02AJF(0.D0) .AND. ABS(CORR2).GT.TOL) THEN
C Since CORR2 did not satisfy the error test we are in a situation
C where CORR2 is `large' and CORR1 is `small'. This indicates
C that Richardson extrapolation is probably not justified.
KXTRAP = JMAX
C WRITE (6,FMT=*) 'Richardson table curtailed (1)'
END IF
IF (CORR2/CORR1.GT.2.D0*TEORAT) THEN
C The last step of Richardson extrapolation was probably
C not justified. Fix KXTRAP accordingly.
KXTRAP = JMAX
C WRITE (6,FMT=*) 'Richardson table curtailed (2)'
END IF
END IF
C WRITE (6,FMT=*) 'KXTRAP=',KXTRAP
IF (KR.GE.MAXTRP) THEN
IFAIL = 8
WRITE (REC,FMT=9290)
9290 FORMAT (' ** More than 10 extrapolations would be required',/,
& ' ** to obtain the required accuracy')
NREC = 2
GO TO 180
END IF
IF (2*NMESH.GT.IWORK) THEN
IFAIL = 11
WRITE (REC,FMT=9300)
9300 FORMAT (' ** IWORK needs to be at least twice as large',/,
& ' ** to obtain the required accuracy')
NREC = 2
GO TO 180
END IF
KR = KR + 1
GO TO 130
C Error handling:
180 IFAIL = P01ABF(IFO,IFAIL,SRNAME,NREC,REC)
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE SL4ARR(A,B,NMESH,XMESH,PP,SP,QP,WP,SL4COF)
C .. Scalar Arguments ..
DOUBLE PRECISION A,B
INTEGER NMESH
C ..
C .. Array Arguments ..
DOUBLE PRECISION PP(NMESH),QP(NMESH),SP(NMESH),WP(NMESH),
& XMESH(0:NMESH)
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4COF
C ..
C .. Local Scalars ..
INTEGER I
C ..
C .. Intrinsic Functions ..
INTRINSIC DBLE
C ..
XMESH(0) = A
DO 10 I = 1,NMESH
XMESH(I) = A + (B-A)*DBLE(I)/DBLE(NMESH)
CALL SL4COF((XMESH(I)+XMESH(I-1))*0.5D0,PP(I),SP(I),QP(I),
& WP(I))
10 CONTINUE
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE SOL4(XMESH,NMESH,IMATCH,P,S,Q,W,ELAM,EPSO,EL,EU,K,TOL,
& KNTMAX,SL4BCS,IFAIL)
C This routine solves a 4th order SL problem with piecewise constant
C coefficients to within a tolerance TOL
C .. Scalar Arguments ..
DOUBLE PRECISION EL,ELAM,EPSO,EU,TOL
INTEGER IFAIL,IMATCH,K,KNTMAX,NMESH
C ..
C .. Array Arguments ..
DOUBLE PRECISION P(NMESH),Q(NMESH),S(NMESH),W(NMESH),
& XMESH(0:NMESH)
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4BCS
C ..
C .. Local Scalars ..
DOUBLE PRECISION DL,DU,E1,E2,EPS,SHIFT,TOLOC
INTEGER IFLAG,IND,KNT
LOGICAL DUAS
C ..
C .. Local Arrays ..
DOUBLE PRECISION C(17)
C ..
C .. External Functions ..
DOUBLE PRECISION SL4MIS
EXTERNAL SL4MIS
C ..
C .. External Subroutines ..
EXTERNAL C05AZF
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,DBLE,MAX,MIN
C ..
IFAIL = 0
EPS = EPSO
SHIFT = DBLE(K) + 5.0D-1
IF (EPS.EQ.0.D0) THEN
EPS = 0.1D0*MAX(1.D0,ABS(ELAM))
END IF
EL = ELAM - EPS
EU = ELAM + EPS
KNT = 0
DUAS = .FALSE.
10 DL = SL4MIS(XMESH,NMESH,IMATCH,P,S,Q,W,EL,SHIFT,SL4BCS,IFAIL)
IF (IFAIL.NE.0) RETURN
C WRITE (6,FMT=*) 'EL, DL:',EL,DL
IF (DL.GT.0.D0) THEN
KNT = KNT + 1
EU = EL
DU = DL
EPS = 2.D0*EPS
EL = EL - EPS
DUAS = .TRUE.
GO TO 10
END IF
IF (DUAS) GO TO 30
DU = SL4MIS(XMESH,NMESH,IMATCH,P,S,Q,W,EU,SHIFT,SL4BCS,IFAIL)
IF (IFAIL.NE.0) RETURN
C WRITE (6,FMT=*) 'EU, DU:',EU,DU
20 IF (DL*DU.GT.0.D0) THEN
C We have not yet bracketed an eigenvalue
KNT = KNT + 1
IF (KNT.GT.KNTMAX) THEN
IFAIL = 9
RETURN
END IF
EPS = 2.D0*EPS
IF (DU.LT.0.D0) THEN
EL = EU
DL = DU
EU = EU + EPS
DU = SL4MIS(XMESH,NMESH,IMATCH,P,S,Q,W,EU,SHIFT,SL4BCS,
& IFAIL)
IF (IFAIL.NE.0) RETURN
C WRITE (6,FMT=*) 'EU, DU (2):',EU,DU
ELSE
EU = EL
DU = DL
EL = EL - EPS
DL = SL4MIS(XMESH,NMESH,IMATCH,P,S,Q,W,EL,SHIFT,SL4BCS,
& IFAIL)
IF (IFAIL.NE.0) RETURN
C WRITE (6,FMT=*) 'EL, DL (2):',EL,DL
END IF
GO TO 20
END IF
C If we are here it means that we have bracketed an eigenvalue and we
C can now think about locating it accurately with C05AZF
30 KNT = 0
IND = -1
TOLOC = TOL
C(1) = DU
IFLAG = 1
E1 = EL
E2 = EU
40 CALL C05AZF(E1,E2,DL,TOLOC,1,C,IND,IFLAG)
IF (IND.EQ.0) GO TO 50
KNT = KNT + 1
IF (KNT.GT.KNTMAX) THEN
IFAIL = 10
RETURN
END IF
DL = SL4MIS(XMESH,NMESH,IMATCH,P,S,Q,W,E1,SHIFT,SL4BCS,IFAIL)
IF (IFAIL.NE.0.D0) RETURN
IF (DL.GT.0.D0) EU = MIN(EU,E1)
IF (DL.LT.0.D0) EL = MAX(EL,E1)
C WRITE (6,FMT=*) 'ELAM D(ELAM),NMESH:',E1,DL,NMESH
IF (ABS(EU-EL).LT.TOLOC*MAX(1.D0,ABS(EL))) GO TO 50
GO TO 40
50 ELAM = E1
RETURN
END
C -------------------------------------------------------------------
DOUBLE PRECISION FUNCTION SL4MIS(XMESH,NMESH,IMATCH,PP,SP,QP,WP,
& ELAM,SHIFT,SL4BCS,IFAIL)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,SHIFT
INTEGER IFAIL,IMATCH,NMESH
C ..
C .. Array Arguments ..
DOUBLE PRECISION PP(1:NMESH),QP(1:NMESH),SP(1:NMESH),WP(1:NMESH),
& XMESH(0:NMESH)
C ..
C .. Local Scalars ..
DOUBLE PRECISION ARGC,OMEGA1,OMEGA2,XEND,XO
INTEGER I,IEND,NL,NR,SIGMA
LOGICAL ISING,PRN
C ..
C .. Local Arrays ..
DOUBLE PRECISION UC(2,2),UL(2,2),UR(2,2),VC(2,2),VL(2,2),VR(2,2)
C ..
C .. External Functions ..
DOUBLE PRECISION X01AAF,X02AJF
EXTERNAL X01AAF,X02AJF
C ..
C .. External Subroutines ..
EXTERNAL CORRCT,F06YAF,SPTH4,ZCOUNT
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4BCS
C ..
C .. Intrinsic Functions ..
INTRINSIC DBLE,MAX,MIN
C ..
IFAIL = 0
PRN = .FALSE.
SL4MIS = 1.D0/X02AJF(0.D0)
C Get the boundary condition at the left-hand endpoint:
CALL SL4BCS(0,ISING,XMESH(0),UL,VL,ELAM)
NL = 0
C Now step across the mesh from XMESH(0) to XMESH(IMATCH),
C carrying the matrices UL, VL such that Theta = (VL+iUL).inv(VL-iUL)
C and carrying ARGL which is Argdet(Theta).
DO 10 I = 1,IMATCH
XO = XMESH(I-1)
XEND = XMESH(I)
IEND = 2
IF (I.EQ.1) IEND = 0
C The routine ZCOUNT does the stepping across the interval
C [XMESH(I-1),XMESH(I)] on which the coefficients P, S, Q and
C W have the constant values PP(I), SP(I), QP(I) and WP(I).
C The eigenparameter has the value ELAM. PRN is a logical
C variable which controls printing of intermediate results
C and IFAIL is an integer which is used to flag errors.
CALL ZCOUNT(UL,VL,PP(I),SP(I),QP(I),WP(I),ELAM,XO,XEND,NL,
& IEND,IFAIL)
IF (IFAIL.NE.0) RETURN
10 CONTINUE
C Get the boundary condition at the right-hand endpoint:
CALL SL4BCS(1,ISING,XMESH(NMESH),UR,VR,ELAM)
NR = 0
DO 20 I = NMESH,IMATCH + 1,-1
IEND = 2
IF (I.EQ.NMESH) IEND = 1
CALL ZCOUNT(UR,VR,PP(I),SP(I),QP(I),WP(I),ELAM,XMESH(I),
& XMESH(I-1),NR,IEND,IFAIL)
IF (IFAIL.NE.0) RETURN
20 CONTINUE
C Compute the central correction sigma:
IEND = 2
IF (IMATCH.EQ.0) IEND = 0
IF (IMATCH.EQ.NMESH) IEND = 1
CALL CORRCT(UL,VL,UR,VR,SIGMA,IEND,IFAIL)
IF (IFAIL.NE.0) THEN
IFAIL = IFAIL + 1
RETURN
END IF
SL4MIS = DBLE(NL+NR+SIGMA) - SHIFT
C Compute the phase angles of Theta_{R}^{*}Theta_{L}
CALL F06YAF('T','N',2,2,2,1.d0,VR,2,UL,2,0.d0,UC,2)
CALL F06YAF('T','N',2,2,2,-1.D0,UR,2,VL,2,1.D0,UC,2)
CALL F06YAF('T','N',2,2,2,1.D0,VR,2,VL,2,0.D0,VC,2)
CALL F06YAF('T','N',2,2,2,1.D0,UR,2,UL,2,1.D0,VC,2)
CALL SPTH4(UC,VC,OMEGA1,OMEGA2,ARGC,'L',IFAIL)
IF (IFAIL.NE.0) THEN
IFAIL = 4
RETURN
END IF
IF (SL4MIS.GE.0.D0 .AND. SL4MIS.LT.1.D0) THEN
SL4MIS = MIN(OMEGA1,OMEGA2)
ELSE IF (SL4MIS.LT.0.D0 .AND. SL4MIS.GT.-1.D0) THEN
SL4MIS = MAX(OMEGA1,OMEGA2)*0.5D0/X01AAF(0.D0) - 1.D0
END IF
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE ZCOUNT(U,V,P,S,Q,W,ELAM,XO,XEND,N,IEND,IFAIL)
C This subroutine integrates the matrices U,V of a 4th order
C Sturm-Liouville problem across an interval on which the
C coefficients are constant, and calculates the nullity count N.
C The formula for the contribution of the current interval to N is
C
C N = NO + SIGMA
C
C where NO is the nullity count across (XO, XEND) of the solution
C UO, satisfying Dirichlet boundary conditions at XEND, and SIGMA
C is the correction parameter CORR at XO, calculated by the
C subroutine CORRCT.
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,P,Q,S,W,XEND,XO
INTEGER IEND,IFAIL,N
C ..
C .. Array Arguments ..
DOUBLE PRECISION U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,B,C,DISC1,DISC2,DUMMY,M,SCAL
INTEGER I,IFLAG,J,NO,SIGMA
LOGICAL LINTSP,LNEG,LPOS,SPCIAL
C ..
C .. Local Arrays ..
DOUBLE PRECISION CU(2,2),CV(2,2),RNORM(2,2),SU(2,2),SV(2,2),
& UO(2,2),UTEMP(2,2),VO(2,2),VTEMP(2,2),Y1(4),
& Y2(4),Y3(4),Y4(4)
C ..
C .. External Subroutines ..
EXTERNAL CORRCT,F06YAF,FUNSOL,INTNLM,LARGLM,NEGLM,SL4CNT
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,MAX,SIN,SQRT
C ..
C Copy the matrices into temporary storage for later
C use.
DO 20 J = 1,2
DO 10 I = 1,2
UTEMP(I,J) = U(I,J)
VTEMP(I,J) = V(I,J)
10 CONTINUE
20 CONTINUE
C
C CALL FUNSOL(XO,XEND,P,S,Q,W,ELAM,Y1,Y2,Y3,Y4,SCAL)
DISC1 = 4.D0*P* (ELAM*W-Q)
DISC2 = S*S + DISC1
LPOS = DISC1 .GE. 0.D0
LINTSP = DISC2 .GE. 0.D0 .AND. S - SQRT(ABS(DISC2)) .GE. 0.D0
LNEG = DISC2 .LT. 0.D0
IF (LPOS .OR. LINTSP) THEN
SCAL = ABS(XEND-XO)*SQRT(0.5D0* (S+SQRT(ABS(DISC2)))/P)
ELSE IF (LNEG) THEN
SCAL = ABS(XEND-XO)*0.5D0*SQRT((S+SQRT(ABS(DISC1)))/P)
ELSE
SCAL = 0.D0
END IF
SPCIAL = SCAL .GT. 1.D-1
IF (.NOT.SPCIAL) THEN
C In this case we must accept the output from FUNSOL as correct
CALL FUNSOL(XO,XEND,P,S,Q,W,ELAM,Y1,Y2,Y3,Y4,SCAL)
CU(1,1) = Y1(1)
CU(1,2) = Y2(1)
CU(2,1) = Y1(2)
CU(2,2) = Y2(2)
SU(1,1) = Y3(1)
SU(1,2) = Y4(1)
SU(2,1) = Y3(2)
SU(2,2) = Y4(2)
SV(1,1) = Y1(3)
SV(1,2) = Y2(3)
SV(2,1) = Y1(4)
SV(2,2) = Y2(4)
CV(1,1) = Y3(3)
CV(1,2) = Y4(3)
CV(2,1) = Y3(4)
CV(2,2) = Y4(4)
C
C Now do the update of U,V from xo to xend, using the formulae
C u(xend) = cu.u(xo)+su.v(xo),
C v(xend) = sv.u(xo)+cv.v(xo).
C
CALL F06YAF('N','N',2,2,2,1.D0,CU,2,UTEMP,2,0.D0,UO,2)
CALL F06YAF('N','N',2,2,2,1.D0,SU,2,VTEMP,2,1.D0,UO,2)
CALL F06YAF('N','N',2,2,2,1.D0,CV,2,VTEMP,2,0.D0,VO,2)
CALL F06YAF('N','N',2,2,2,1.D0,SV,2,UTEMP,2,1.D0,VO,2)
C
C Now do some postmultiplication to avoid ill-conditioning
C of these matrices.
C
B = -UO(1,1)*VO(2,2) + UO(1,2)*VO(2,1) - UO(2,2)*VO(1,1) +
& UO(2,1)*VO(1,2)
AC = UO(1,1)*UO(2,2) - UO(2,1)*UO(1,2) - VO(1,1)*VO(2,2) +
& VO(2,1)*VO(1,2)
C = UO(1,1)*UO(2,2) - UO(2,1)*UO(1,2) + VO(1,1)*VO(2,2) -
& VO(2,1)*VO(1,2)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C = COS(M/2.D0)
AC = SIN(M/2.D0)
C Form the normalising matrix:
DO 40 J = 1,2
DO 30 I = 1,2
RNORM(I,J) = UO(I,J)*C - VO(I,J)*AC
30 CONTINUE
40 CONTINUE
C Form the inverse of the normalising matrix:
C = RNORM(1,1)*RNORM(2,2) - RNORM(1,2)*RNORM(2,1)
AC = RNORM(1,1)/C
RNORM(1,1) = RNORM(2,2)/C
RNORM(2,2) = AC
RNORM(1,2) = -RNORM(1,2)/C
RNORM(2,1) = -RNORM(2,1)/C
C Peform the postmultiplication: U = UO.RNORM, V = VO.RNORM
CALL F06YAF('N','N',2,2,2,1.D0,UO,2,RNORM,2,0.D0,U,2)
CALL F06YAF('N','N',2,2,2,1.D0,VO,2,RNORM,2,0.D0,V,2)
C
C Now set the U and V (evaluated at X = 0) for Dirichlet conditions
C at X = XEND. (Thus U(XEND) = 0, V(XEND) = I.)
UO(1,1) = -SU(1,1)
UO(1,2) = SU(1,2)
UO(2,1) = SU(2,1)
UO(2,2) = -SU(2,2)
VO(1,1) = CV(1,1)
VO(1,2) = -CV(1,2)
VO(2,1) = -CV(2,1)
VO(2,2) = CV(2,2)
ELSE
IF (LPOS) THEN
CALL LARGLM(ELAM,P,S,Q,W,U,V,UO,VO,XO,XEND)
ELSE IF (LINTSP) THEN
CALL INTNLM(ELAM,P,S,Q,W,U,V,UO,VO,XO,XEND)
ELSE IF (LNEG) THEN
CALL NEGLM(ELAM,P,S,Q,W,U,V,UO,VO,XO,XEND)
END IF
END IF
C Normalise U and V so that the greatest element of the two is
C equal to 1.
DUMMY = 0.D0
DO 60 J = 1,2
DO 50 I = 1,2
DUMMY = MAX(DUMMY,ABS(U(I,J)),ABS(V(I,J)))
50 CONTINUE
60 CONTINUE
DO 80 J = 1,2
DO 70 I = 1,2
U(I,J) = U(I,J)/DUMMY
V(I,J) = V(I,J)/DUMMY
70 CONTINUE
80 CONTINUE
C
C We have now updated U, V from XO to XEND.
C Now we will calculate the correction parameter SIGMA
C
IF (XO.LE.XEND) THEN
CALL CORRCT(UTEMP,VTEMP,UO,VO,SIGMA,IEND,IFLAG)
ELSE
CALL CORRCT(UO,VO,UTEMP,VTEMP,SIGMA,IEND,IFLAG)
END IF
IF (IFLAG.NE.0) THEN
IF (IFLAG.EQ.1) IFAIL = 6
IF (IFLAG.EQ.2) IFAIL = 7
RETURN
END IF
C
C Next find the nullity count for UO.
C
CALL SL4CNT(XEND-XO,P,S,Q,W,ELAM,NO)
C
C Now we can update the total nullity count N:
C
N = N + NO + SIGMA
END
C --------------------------------------------------------------------
SUBROUTINE CORRCT(UL,VL,UR,VR,CORR,IEND,IFLAG)
C
C This subroutine calculates the correction term "corr"
C in the formula N(lambda) = No(lambda) + corr.
C The correction term is calculated at an endpoint xo of
C a mesh interval. At this point we have input data (UL,VL)
C and (UR,VR). There are 4 cases to consider, depending on
C whether det(UL) and det(UR) are or are not zero. It is convenient
C to consider a 5th case, when UL or UR is the zero matrix. Thus,
C we have the following cases:
C (1) det(UL) and det(UR) are not zero.
C (2) One of the matrices UL or UR is zero.
C (3) det(UL) = 0; UL and det(UR) are not zero.
C (4) det(UL) and UR are not zero; det(UR) = 0.
C (5) det(UL) = det(UR) =0; UL and UR are not zero.
C The integer IFLAG is normally 0, but equals 1 if
C det(DIFF) > 0, while trace(DIFF) = 0,
C where DIFF = WL - WR, and equals 2 if the matrices U
C and V do not fit any of the cases.
C
C IEND tells CORRCT if it is being called at the end of an
C interval.
C IEND = 0 means left endpoint
C IEND = 1 means right endpoint
C IEND = 2 means not an endpoint.
C
C--------------------------------------------------------------------
C
C
C .. Scalar Arguments ..
INTEGER CORR,IEND,IFLAG
C ..
C .. Array Arguments ..
DOUBLE PRECISION UL(2,2),UR(2,2),VL(2,2),VR(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION A,ADL,ADR,AL,AR,DET,DUL,DUR,DVL,DVR,F,G,KL,KR,ML,
& MR,NL,NR,SGN,SGNL,SGNR,SL,SR,TRC
LOGICAL KN,MN,ZL,ZR
C ..
C .. Local Arrays ..
DOUBLE PRECISION DIFF(2,2),VUAL(2,2),VUAR(2,2)
C ..
C .. External Subroutines ..
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,SIGN
C ..
IFLAG = 0
DUL = UL(1,1)*UL(2,2) - UL(1,2)*UL(2,1)
DUR = UR(1,1)*UR(2,2) - UR(1,2)*UR(2,1)
C
DVL = VL(1,1)*VL(2,2) - VL(1,2)*VL(2,1)
DVR = VR(1,1)*VR(2,2) - VR(1,2)*VR(2,1)
C
VUAL(1,1) = UL(2,2)*VL(1,1) - UL(2,1)*VL(1,2)
VUAL(1,2) = UL(2,2)*VL(2,1) - UL(2,1)*VL(2,2)
VUAL(2,1) = VUAL(1,2)
VUAL(2,2) = UL(1,1)*VL(2,2) - UL(1,2)*VL(2,1)
C
VUAR(1,1) = UR(2,2)*VR(1,1) - UR(2,1)*VR(1,2)
VUAR(1,2) = UR(2,2)*VR(2,1) - UR(2,1)*VR(2,2)
VUAR(2,1) = VUAR(1,2)
VUAR(2,2) = UR(1,1)*VR(2,2) - UR(1,2)*VR(2,1)
C
SGNL = SIGN(1.D0,DUL)
SGNR = SIGN(1.D0,DUR)
C
C Case (1): det(UL) and det(UR) are not zero.
IF (DUL.NE.0.D0 .AND. DUR.NE.0.D0) THEN
DIFF(1,1) = VUAL(1,1)*ABS(DUR)*SGNL - VUAR(1,1)*ABS(DUL)*SGNR
DIFF(1,2) = VUAL(1,2)*ABS(DUR)*SGNL - VUAR(1,2)*ABS(DUL)*SGNR
DIFF(2,1) = DIFF(1,2)
DIFF(2,2) = VUAL(2,2)*ABS(DUR)*SGNL - VUAR(2,2)*ABS(DUL)*SGNR
C
DET = DIFF(1,1)*DIFF(2,2) - DIFF(1,2)*DIFF(2,1)
TRC = DIFF(1,1) + DIFF(2,2)
IF (DET.LT.0.D0) THEN
CORR = 1
RETURN
C
ELSE IF (DET.GT.0.D0) THEN
IF (TRC.LT.0.D0) THEN
CORR = 2
RETURN
ELSE IF (TRC.GT.0.D0) THEN
CORR = 0
RETURN
ELSE IF (TRC.EQ.0.D0) THEN
IFLAG = 1
RETURN
END IF
ELSE IF (DET.EQ.0.D0) THEN
IF (TRC.LT.0.D0) THEN
CORR = 1
RETURN
ELSE
CORR = 0
RETURN
END IF
END IF
END IF
C
C
C Case (2): UL = 0 or UR = 0.
C
C
ZL = UL(1,1) .EQ. 0.D0 .AND. UL(1,2) .EQ. 0.D0 .AND.
& UL(2,1) .EQ. 0.D0 .AND. UL(2,2) .EQ. 0.D0
ZR = UR(1,1) .EQ. 0.D0 .AND. UR(1,2) .EQ. 0.D0 .AND.
& UR(2,1) .EQ. 0.D0 .AND. UR(2,2) .EQ. 0.D0
C
C
IF (ZL .OR. ZR) THEN
IF (IEND.EQ.0) THEN
IF (ZL) THEN
CORR = 0
RETURN
ELSE IF (DUL.EQ.0.D0) THEN
CORR = 1
RETURN
ELSE
CORR = 2
RETURN
END IF
ELSE IF (IEND.EQ.1) THEN
IF (ZR) THEN
CORR = 0
RETURN
ELSE IF (DUR.EQ.0.D0) THEN
CORR = 1
RETURN
ELSE
CORR = 2
RETURN
END IF
ELSE IF (IEND.EQ.2) THEN
CORR = 2
RETURN
END IF
END IF
C
C
KL = VUAL(1,1)
ML = VUAL(1,2)
NL = VUAL(2,2)
C
KR = VUAR(1,1)
MR = VUAR(1,2)
NR = VUAR(2,2)
C
C
C Case (3): det(UL) = 0; UL and det(UR) are not zero.
C
C
IF (DUL.EQ.0 .AND. DUR.NE.0) THEN
SGN = SIGN(1.D0,KL+NL)
A = ABS(KL+NL)
ADR = ABS(DUR)
IF (KL.NE.0.D0) THEN
F = (KR*ML*ML-2*MR*KL*ML+NR*KL*KL)*A*SGNR
G = DVL* (KL*KL+ML*ML)*ADR*SGN
IF (F.LE.G) THEN
IF (IEND.EQ.0) THEN
CORR = 0
RETURN
ELSE
CORR = 1
RETURN
END IF
ELSE
IF (IEND.EQ.0) THEN
CORR = 1
RETURN
ELSE
CORR = 2
RETURN
END IF
END IF
ELSE
F = KR*A*SGNR
G = DVL*ADR*SGN
IF (F.LE.G) THEN
IF (IEND.EQ.0) THEN
CORR = 0
RETURN
ELSE
CORR = 1
RETURN
END IF
ELSE
IF (IEND.EQ.0) THEN
CORR = 1
RETURN
ELSE
CORR = 2
RETURN
END IF
END IF
END IF
END IF
C
C
C Case (4): det(UL) and UR are not zero; det (UR) = 0.
C
C
IF (DUL.NE.0 .AND. DUR.EQ.0) THEN
SGN = SIGN(1.D0,KR+NR)
A = ABS(KR+NR)
ADL = ABS(DUL)
IF (KR.NE.0.D0) THEN
F = (KL*MR*MR-2*ML*KR*MR+NL*KR*KR)*A*SGNL
G = DVR* (KR*KR+MR*MR)*ADL*SGN
IF (F.GE.G) THEN
IF (IEND.EQ.1) THEN
CORR = 0
RETURN
ELSE
CORR = 1
RETURN
END IF
ELSE
IF (IEND.EQ.1) THEN
CORR = 1
RETURN
ELSE
CORR = 2
RETURN
END IF
END IF
ELSE
F = KL*A*SGNL
G = DVR*ADL*SGN
IF (F.GE.G) THEN
IF (IEND.EQ.1) THEN
CORR = 0
RETURN
ELSE
CORR = 1
RETURN
END IF
ELSE
IF (IEND.EQ.1) THEN
CORR = 1
RETURN
ELSE
CORR = 2
RETURN
END IF
END IF
END IF
END IF
C
C
C Case (5): det(UL) = det(UR) = 0; UL and UR are not zero.
C
C
IF (DUL.EQ.0.D0 .AND. DUR.EQ.0.D0) THEN
KN = KL*NR .EQ. KR*NL
MN = ML*NR .EQ. MR*NL
AL = ABS(KL+NL)
SL = SIGN(1.D0,KL+NL)
AR = ABS(KR+NR)
SR = SIGN(1.D0,KR+NR)
F = DVL*AR*SL
G = DVR*AL*SR
IF (KN .AND. MN .AND. F.GE.G) THEN
IF (IEND.EQ.2) THEN
CORR = 1
RETURN
ELSE
CORR = 0
RETURN
END IF
ELSE
IF (IEND.EQ.2) THEN
CORR = 2
RETURN
ELSE
CORR = 1
RETURN
END IF
END IF
END IF
IFLAG = 2
RETURN
END
C --------------------------------------------------------------------
SUBROUTINE LARGLM(ELAM,P,S,Q,W,U,V,SU,CV,XO,XEND)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,P,Q,S,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION CV(2,2),SU(2,2),U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,AOMX,B,BETA,C,COSHCS,COSHNX,COSHSN,COSOMX,
& DETCU,DETCV,DETSU,DETSV,DETU,DETV,DISC,FAC,H,H2,
& H2SNSN,M,NO,NO2,NU,NU2,NU3,NU4,NU6,NUX,OMEGA,
& OMEGA2,OMEGA3,OMEGA4,OMEGA6,OMX,PNO,PNO2,PNORAD,
& PRAD,PRAD2,RAD,RAD2,SINHCS,SINHNX,SINHSN,SINOMX,
& TRUVA
INTEGER I,J
C ..
C .. Local Arrays ..
DOUBLE PRECISION A(2,2),CUSVA(2,2),SUCVA(2,2),T(2,2),UUA(2,2),
& UVA(2,2),VUA(2,2),VVA(2,2)
C ..
C .. External Functions ..
DOUBLE PRECISION PHI
EXTERNAL PHI
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,EXP,SIN,SQRT
C ..
C .. External Subroutines ..
C ..
DISC = SQRT(S*S+4.D0*P* (ELAM*W-Q))
NU = SQRT(0.5d0* (S+DISC)/P)
OMEGA = SQRT(0.5d0* (-S+DISC)/P)
DETU = U(1,1)*U(2,2) - U(2,1)*U(1,2)
DETV = V(1,1)*V(2,2) - V(2,1)*V(1,2)
H = XEND - XO
H2 = H*H
NUX = ABS(NU*H)
FAC = EXP(-NUX)
IF (NUX.GE.0.1D0) THEN
SINHNX = (1.D0-FAC*FAC)/ (2.D0*NUX)
ELSE
SINHNX = PHI(NUX*NUX)*FAC
END IF
COSHNX = (1.D0+FAC*FAC)/2.D0
OMX = OMEGA*H
AOMX = ABS(OMX)
IF (AOMX.GE.0.1D0) THEN
SINOMX = SIN(OMX)/OMX
ELSE
SINOMX = PHI(-OMX*OMX)
END IF
COSOMX = COS(OMX)
COSHSN = H*COSHNX*SINOMX
SINHCS = H*SINHNX*COSOMX
COSHCS = COSHNX*COSOMX
SINHSN = SINHNX*SINOMX
H2SNSN = H2*SINHSN
OMEGA2 = OMEGA*OMEGA
NU2 = NU*NU
OMEGA3 = OMEGA*OMEGA2
NU3 = NU*NU2
OMEGA4 = OMEGA2*OMEGA2
NU4 = NU2*NU2
OMEGA6 = OMEGA2*OMEGA4
NU6 = NU2*NU4
RAD = OMEGA2 + NU2
RAD2 = RAD*RAD
PRAD = P*RAD
PRAD2 = PRAD*PRAD
NO = NU*OMEGA
PNO = P*NO
PNORAD = PNO*RAD
NO2 = NO*NO
PNO2 = P*NO2
C Compute SU.Adjoint(CV).exp(-NUX):
SUCVA(1,1) = (COSHSN-SINHCS)/PRAD
SUCVA(1,2) = ((OMEGA2-NU2)* (FAC-COSHCS)+2.D0*NO2*H2SNSN)/
& (PRAD*RAD)
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = (NU2*SINHCS+OMEGA2*COSHSN)/PRAD
C Compute CU.Adjoint(SV).exp(-NUX):
CUSVA(1,1) = P* (NU4*SINHCS-OMEGA4*COSHSN)/RAD
CUSVA(1,2) = -PNO2* ((NU4+OMEGA4)*H2SNSN+
& (NU2-OMEGA2)* (FAC-COSHCS))/RAD2
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = -PNO2* (NU2*COSHSN+OMEGA2*SINHCS)/RAD
C Compute A = U.Adjoint(V):
A(1,1) = U(1,1)*V(2,2) - U(1,2)*V(2,1)
A(1,2) = U(1,2)*V(1,1) - U(1,1)*V(1,2)
A(2,1) = A(1,2)
A(2,2) = V(1,1)*U(2,2) - U(2,1)*V(1,2)
C Compute T = (CU.A.Adjoint(CV) + SU.A.Adjoint(SV))*exp(-NUX):
T(1,1) = A(1,1)*COSHCS + 2.D0*A(1,2)* (OMEGA2*COSHSN+NU2*SINHCS)/
& RAD + A(2,2)*H2SNSN
T(1,2) = (A(1,1)*NO2* (SINHCS-COSHSN)+
& A(1,2)* ((OMEGA2-NU2)* (FAC* (OMEGA2-NU2)+
& 2.D0*NO2*H2SNSN)+4.D0*NO2*COSHCS)/RAD+
& A(2,2)* (NU2*SINHCS+OMEGA2*COSHSN))/RAD
T(2,1) = T(1,2)
T(2,2) = -A(1,1)*NO2*H2SNSN + 2.D0*A(1,2)*NO2* (SINHCS-COSHSN)/
& RAD + A(2,2)*COSHCS
C Compute U.Adjoint(V) = (CU.Uo+SU.Vo).Adjoint(SV.Uo+CV.Vo):
UVA(1,1) = DETU*CUSVA(1,1) + DETV*SUCVA(1,1) + T(1,1)
UVA(1,2) = DETU*CUSVA(1,2) + DETV*SUCVA(1,2) + T(1,2)
UVA(2,1) = DETU*CUSVA(2,1) + DETV*SUCVA(2,1) + T(2,1)
UVA(2,2) = DETU*CUSVA(2,2) + DETV*SUCVA(2,2) + T(2,2)
VUA(1,1) = UVA(2,2)
VUA(2,2) = UVA(1,1)
VUA(1,2) = -UVA(1,2)
VUA(2,1) = -UVA(2,1)
DETCU = ((NU4+OMEGA4)*COSHCS+NO2* ((NU2-OMEGA2)*H2SNSN+2.D0*FAC))/
& RAD2
DETSU = ((NU2-OMEGA2)*H2SNSN+2.D0* (FAC-COSHCS))/PRAD2
DETCV = DETCU
DETSV = P*PNO* ((OMEGA6-NU6)*NO*H2SNSN+
& 2.D0*NU3*OMEGA3* (FAC-COSHCS))/RAD2
C Compute trace(CU.A.Adjoint(SU)).exp(-NUX):
BETA = A(1,1)* (OMEGA2*COSHSN+NU2*SINHCS)/PRAD -
& 2.D0*A(1,2)* ((NU2-OMEGA2)* (FAC-COSHCS)-2.D0*NO2*H2SNSN)/
& (P*RAD2) + A(2,2)* (COSHSN-SINHCS)/PRAD
UUA(1,1) = DETU*DETCU + DETV*DETSU + BETA
UUA(1,2) = 0.D0
UUA(2,1) = 0.D0
UUA(2,2) = UUA(1,1)
C Compute trace(CV.A.Adjoint(SV)).exp(-NUX):
BETA = (-A(1,1)*PNO2* (NU2*COSHSN+OMEGA2*SINHCS)-
& 2.D0*PNO2*A(1,2)* ((NU2-OMEGA2)* (FAC-COSHCS)+ (NU4+
& OMEGA4)*H2SNSN)/RAD-A(2,2)*P* (OMEGA4*COSHSN-NU4*SINHCS))/
& RAD
VVA(1,1) = DETU*DETSV + DETV*DETCV + BETA
VVA(1,2) = 0.D0
VVA(2,1) = 0.D0
VVA(2,2) = VVA(1,1)
TRUVA = UVA(1,1) + UVA(2,2)
AC = UUA(1,1) - VVA(1,1)
B = -TRUVA
C = UUA(1,1) + VVA(1,1)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
DO 20 J = 1,2
DO 10 I = 1,2
U(I,J) = UUA(I,J)*COS(M/2.D0) - UVA(I,J)*SIN(M/2.D0)
V(I,J) = VUA(I,J)*COS(M/2.D0) - VVA(I,J)*SIN(M/2.D0)
10 CONTINUE
20 CONTINUE
C That's done the very tricky renormalisation. Now do the SU and CV
C matrices which are the U and V matrices for Theta_{0}.
TRUVA = SUCVA(1,1) + SUCVA(2,2)
AC = (DETSU-DETCV)
B = -TRUVA
C = (DETSU+DETCV)
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C NOTE: in the following formulae we want CV and SU for BACKWARDS
C integration. We exploit the fact that the diagonal terms of
C SUCVA are odd functions of the direction of integration while
C the off-diagonal terms are even functions of the direction
C of integration, while the determinants which appear in the
C formulae are even functions of the direction of integration.
SU(1,1) = DETSU*COS(M/2.D0) + SUCVA(1,1)*SIN(M/2.D0)
SU(1,2) = -SUCVA(1,2)*SIN(M/2.D0)
SU(2,1) = -SUCVA(2,1)*SIN(M/2.D0)
SU(2,2) = DETSU*COS(M/2.D0) + SUCVA(2,2)*SIN(M/2.D0)
CV(1,1) = -SUCVA(2,2)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
CV(1,2) = -SUCVA(1,2)*COS(M/2.D0)
CV(2,1) = -SUCVA(2,1)*COS(M/2.D0)
CV(2,2) = -SUCVA(1,1)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE INTNLM(ELAM,P,S,Q,W,U,V,SU,CV,XO,XEND)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,P,Q,S,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION CV(2,2),SU(2,2),U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,AEHH,AETAH,ARHH,B,BETA,C,CMSN,CNCM,CNSM,
& COSHMX,COSHNX,DETCU,DETCV,DETSU,DETSV,DETU,DETV,
& DISC,ETA,ETA2,ETAH,ETASQ,FAC,FACL,H,H2,M,MU,MU2,
& MU3,MU4,MU6,MUX,NM,NU,NU2,NU3,NU4,NU6,NUX,PNM,
& PNMR,PRAD,PRAD2,RAD,RAD2,RADH2,RH,RH2,RHH,RHSQ,
& SINHMX,SINHNX,SN2ERT,SN2HRT,SNERAT,SNHR1,SNHR2,
& SNHRAT,SNSM,TRUVA
INTEGER I,J
LOGICAL SMALL
C ..
C .. Local Arrays ..
DOUBLE PRECISION A(2,2),CUSVA(2,2),SUCVA(2,2),T(2,2),UUA(2,2),
& UVA(2,2),VUA(2,2),VVA(2,2)
C ..
C .. External Functions ..
DOUBLE PRECISION PHI,SNHDIF
EXTERNAL PHI,SNHDIF
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,EXP,MIN,SIGN,SIN,SQRT
C ..
DISC = SQRT(S*S+4.D0*P* (ELAM*W-Q))
NU2 = 0.5D0* (S+DISC)/P
MU2 = NU2 - DISC/P
NU = SQRT(NU2)
MU = SQRT(MU2)
NU3 = NU*NU2
MU3 = MU*MU2
NU4 = NU2*NU2
MU4 = MU2*MU2
NU6 = NU3*NU3
MU6 = MU3*MU3
DETU = U(1,1)*U(2,2) - U(2,1)*U(1,2)
DETV = V(1,1)*V(2,2) - V(2,1)*V(1,2)
H = XEND - XO
H2 = H*H
NUX = ABS(NU*H)
SINHNX = (1.D0-EXP(-2.D0*NUX))/2.D0
COSHNX = 1.D0 - SINHNX
SINHNX = SINHNX*SIGN(1.D0,H)
MUX = ABS(MU*H)
SINHMX = (1.D0-EXP(-2.D0*MUX))/2.D0
COSHMX = 1.D0 - SINHMX
SINHMX = SINHMX*SIGN(1.D0,H)
FAC = EXP(-NUX-MUX)
CMSN = COSHMX*SINHNX
CNSM = COSHNX*SINHMX
CNCM = COSHNX*COSHMX
SNSM = SINHNX*SINHMX
C If any of the following quantities is small then we are
C in trouble; Taylor series expansions have to be used.
C RAD = NU**2-MU**2
RAD = DISC/P
RAD2 = RAD*RAD
PRAD = DISC
PRAD2 = DISC*DISC
NM = NU*MU
PNM = P*NM
PNMR = PNM*RAD
SMALL = H2*MIN(ABS(RAD),ABS(PRAD),ABS(PNMR),PNM) .LT. 1.D-6
C Compute SU.Adjoint(CV).exp(-(NUX+MUX)):
IF (SMALL) THEN
C Compute certain quantities which will be jolly useful
C when calculating the matrix elements
ETA2 = NU + MU
ETA = 0.5D0*ETA2
ETASQ = ETA*ETA
ETAH = ETA*H
AETAH = ABS(ETAH)
RH2 = NU - MU
RH = 0.5D0*RH2
RHSQ = RH*RH
RHH = RH*H
ARHH = ABS(RHH)
IF (ABS(RHH).LT.0.1D0) THEN
SNHRAT = PHI(RHH*RHH)
SNHRAT = SNHRAT*SNHRAT*FAC
SN2HRT = PHI(4.D0*RHH*RHH)*FAC
ELSE
FACL = EXP(-MUX)
SNHRAT = 0.5D0* (1.D0-EXP(-2.D0*ARHH))/ARHH
SNHRAT = SNHRAT*FACL
SNHRAT = SNHRAT*SNHRAT
SN2HRT = (1.D0-EXP(-4.D0*ARHH))/ (4.D0*ARHH)
SN2HRT = SN2HRT*FACL*FACL
END IF
IF (ABS(ETAH).LT.0.1D0) THEN
SNERAT = PHI(ETAH*ETAH)
SNERAT = SNERAT*SNERAT*FAC
SN2ERT = PHI(4.D0*ETAH*ETAH)*FAC
ELSE
SNERAT = 0.5D0* (1.D0-EXP(-2.D0*AETAH))/AETAH
SNERAT = SNERAT*SNERAT
SN2ERT = (1.D0-EXP(-4.D0*AETAH))/ (4.D0*AETAH)
END IF
IF (MUX.LT.0.1D0) THEN
SNHR1 = PHI(MUX*MUX)*EXP(-MUX)
ELSE
SNHR1 = SIGN(1.D0,H)*SINHMX/MUX
END IF
IF (NUX.LT.0.1D0) THEN
SNHR2 = PHI(NUX*NUX)*EXP(-NUX)
ELSE
SNHR2 = SIGN(1.D0,H)*SINHNX/NUX
END IF
END IF
IF (.NOT.SMALL) THEN
SUCVA(1,1) = (NU*CNSM-MU*CMSN)/PNMR
SUCVA(1,2) = - (S* (FAC-CNCM)+2.D0*PNM*SNSM)/PRAD2
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = (NU*CMSN-MU*CNSM)/PRAD
ELSE
RADH2 = ABS(RAD)*H2
AEHH = ABS(ETASQ-RHSQ)*H2
IF (RADH2.GE.0.01D0 .AND. AEHH.LT.0.01D0) THEN
SUCVA(1,1) = H* (COSHNX*SNHR1-COSHMX*SNHR2)/PRAD
ELSE IF (RADH2.LT.0.01D0 .AND. AEHH.GE.0.01D0) THEN
SUCVA(1,1) = H* (SN2ERT-SN2HRT)/ (2.D0*PNM)
ELSE
SUCVA(1,1) = 2.D0* (H**3)*FAC*
& SNHDIF(2.D0*AETAH,2.D0*ARHH)/P
END IF
SUCVA(1,2) = 0.25D0*H*H* (SNERAT+SNHRAT)/P
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = 0.5D0*H* (SN2ERT+SN2HRT)/P
END IF
C Compute CU.Adjoint(SV).exp(-NUX-MUX):
IF (.NOT.SMALL) THEN
CUSVA(1,1) = P* (NU3*CMSN-MU3*CNSM)/RAD
CUSVA(1,2) = PNM* (NM* (NU2+MU2)* (FAC-CNCM)+ (NU4+MU4)*SNSM)/
& RAD2
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = PNM* (NU3*CNSM-MU3*CMSN)/RAD
ELSE
FACL = 0.25D0*PNM*H
CUSVA(1,1) = P*H*0.5D0* ((ETASQ+3.D0*RHSQ)*SN2HRT+
& (RHSQ+3.D0*ETASQ)*SN2ERT)
CUSVA(1,2) = 0.25D0*PNM* (H2* (RHSQ*SNERAT-ETASQ*SNHRAT)+
& 3.D0*SNSM)
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = 2.D0*FACL* ((RHSQ+3.D0*ETASQ)*SN2ERT-
& (ETASQ+3.D0*RHSQ)*SN2HRT)
END IF
C Compute A = U.Adjoint(V):
A(1,1) = U(1,1)*V(2,2) - U(1,2)*V(2,1)
A(1,2) = U(1,2)*V(1,1) - U(1,1)*V(1,2)
A(2,1) = A(1,2)
A(2,2) = V(1,1)*U(2,2) - U(2,1)*V(1,2)
C Compute T = (CU.A.Adjoint(CV) + SU.A.Adjoint(SV))*exp(-NUX-MUX):
IF (.NOT.SMALL) THEN
T(1,1) = A(1,1)*CNCM + 2.D0*A(1,2)* (NU*CMSN-MU*CNSM)/RAD +
& A(2,2)*SNSM/NM
T(1,2) = (A(1,1)*NM* (NU*CNSM-MU*CMSN)+
& A(1,2)* ((FAC* (NU2+MU2)+2.D0*NM*SNSM)* (NU2+MU2)-
& 4.D0*NM*NM*CNCM)/RAD+A(2,2)* (NU*CMSN-MU*CNSM))/RAD
T(2,1) = T(1,2)
T(2,2) = NM* (A(1,1)*SNSM+2.D0*A(1,2)* (NU*CNSM-MU*CMSN)/
& RAD) + A(2,2)*CNCM
ELSE
T(1,1) = A(1,1)*CNCM + 2.D0*P*A(1,2)*SUCVA(2,2) +
& A(2,2)*H2*SNHR1*SNHR2
T(1,2) = A(1,1)*PNM*NM*SUCVA(1,1) +
& 0.5D0*A(1,2)* (FAC+CNCM-H2*
& (ETASQ*SNHRAT+RHSQ*SNERAT)) + A(2,2)*P*SUCVA(2,2)
T(2,1) = T(1,2)
T(2,2) = NM* (A(1,1)*SNSM+2.D0*A(1,2)*PNM*SUCVA(1,1)) +
& A(2,2)*CNCM
END IF
C Compute U.Adjoint(V) = (CU.Uo+SU.Vo).Adjoint(SV.Uo+CV.Vo):
UVA(1,1) = DETU*CUSVA(1,1) + DETV*SUCVA(1,1) + T(1,1)
UVA(1,2) = DETU*CUSVA(1,2) + DETV*SUCVA(1,2) + T(1,2)
UVA(2,1) = DETU*CUSVA(2,1) + DETV*SUCVA(2,1) + T(2,1)
UVA(2,2) = DETU*CUSVA(2,2) + DETV*SUCVA(2,2) + T(2,2)
VUA(1,1) = UVA(2,2)
VUA(2,2) = UVA(1,1)
VUA(1,2) = -UVA(1,2)
VUA(2,1) = -UVA(2,1)
IF (.NOT.SMALL) THEN
DETCU = ((NU4+MU4)*CNCM-NM* ((NU2+MU2)*SNSM+2.D0*NM*FAC))/RAD2
DETCV = DETCU
DETSU = (2.D0*NM* (FAC-CNCM)+ (NU2+MU2)*SNSM)/ (PRAD2*NM)
DETSV = P*P*NM* (2.D0*NU3*MU3* (FAC-CNCM)+ (NU6+MU6)*SNSM)/
& RAD2
ELSE
DETCU = 0.25D0* (H2* (ETASQ*SNHRAT+RHSQ*SNERAT)+FAC+3.D0*CNCM)
DETCV = DETCU
RADH2 = ABS(RAD)*H2
AEHH = ABS(ETASQ-RHSQ)*H2
IF (RADH2.GE.0.01D0 .AND. AEHH.LT.0.01D0) THEN
FACL = H2*SNHR1*SNHR2
DETSU = (2.D0* (FAC-CNCM)+ (NU2+MU2)*FACL)/PRAD2
ELSE IF (RADH2.LT.0.01D0 .AND. AEHH.GE.0.01D0) THEN
FACL = H*0.5D0/P
DETSU = (SNERAT-SNHRAT)*FACL*FACL/ (ETASQ-RHSQ)
ELSE
FACL = 0.5D0*H2/P
DETSU = (SQRT(SNERAT)-SQRT(SNHRAT))*SNHDIF(AETAH,ARHH)*
& FACL*FACL
END IF
DETSV = (P*PNM/8.D0)* (2.D0*H*H*
& (RHSQ*RHSQ*SNERAT-ETASQ*ETASQ*SNHRAT)+
& 3.D0*NM* (CNCM-FAC)+15.D0* (ETASQ+RHSQ)*SNSM)
END IF
C Compute trace(CU.A.Adjoint(SU)).exp(-NUX-MUX):
C IF (.NOT.SMALL) THEN
C
C BETA = A(1,1)* (NU*CMSN-MU*CNSM)/PRAD +
C & 2.D0*A(1,2)* ((NU2+MU2)* (CNCM-FAC)-2.D0*NM*SNSM)/
C & (P*RAD2) + A(2,2)* (NU*CNSM-MU*CMSN)/PNMR
C
C ELSE
BETA = A(1,1)*SUCVA(2,2) + 2.D0*A(1,2)*SUCVA(1,2) +
& A(2,2)*SUCVA(1,1)
C END IF
UUA(1,1) = DETU*DETCU + DETV*DETSU + BETA
UUA(1,2) = 0.D0
UUA(2,1) = 0.D0
UUA(2,2) = UUA(1,1)
C Compute trace(CV.A.Adjoint(SV)).exp(-NUX-MUX):
C IF (.NOT.SMALL) THEN
C
C BETA = (A(1,1)*PNM* (NU3*CNSM-MU3*CMSN)+
C & 2.D0*A(1,2)*PNM* ((NU4+MU4)*SNSM+NM* (NU2+MU2)* (FAC-
C & CNCM))/RAD+A(2,2)*P* (NU3*CMSN-MU3*CNSM))/RAD
C
C ELSE
BETA = A(1,1)*CUSVA(2,2) + 2.D0*A(1,2)*CUSVA(1,2) +
& A(2,2)*CUSVA(1,1)
C END IF
VVA(1,1) = DETU*DETSV + DETV*DETCV + BETA
VVA(1,2) = 0.D0
VVA(2,1) = 0.D0
VVA(2,2) = VVA(1,1)
TRUVA = UVA(1,1) + UVA(2,2)
AC = UUA(1,1) - VVA(1,1)
B = -TRUVA
C = UUA(1,1) + VVA(1,1)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
DO 20 J = 1,2
DO 10 I = 1,2
U(I,J) = UUA(I,J)*COS(M/2.D0) - UVA(I,J)*SIN(M/2.D0)
V(I,J) = VUA(I,J)*COS(M/2.D0) - VVA(I,J)*SIN(M/2.D0)
10 CONTINUE
20 CONTINUE
C That's done the very tricky renormalisation. Now do the SU and CV
C matrices which are the U and V matrices for Theta_{0}.
TRUVA = SUCVA(1,1) + SUCVA(2,2)
AC = (DETSU-DETCV)
B = -TRUVA
C = (DETSU+DETCV)
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C NOTE: in the following formulae we want CV and SU for BACKWARDS
C integration. We exploit the fact that the diagonal terms of
C SUCVA are odd functions of the direction of integration while
C the off-diagonal terms are even functions of the direction
C of integration, while the determinants which appear in the
C formulae are even functions of the direction of integration.
SU(1,1) = DETSU*COS(M/2.D0) + SUCVA(1,1)*SIN(M/2.D0)
SU(1,2) = -SUCVA(1,2)*SIN(M/2.D0)
SU(2,1) = -SUCVA(2,1)*SIN(M/2.D0)
SU(2,2) = DETSU*COS(M/2.D0) + SUCVA(2,2)*SIN(M/2.D0)
CV(1,1) = -SUCVA(2,2)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
CV(1,2) = -SUCVA(1,2)*COS(M/2.D0)
CV(2,1) = -SUCVA(2,1)*COS(M/2.D0)
CV(2,2) = -SUCVA(1,1)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE NEGLM(ELAM,P,S,Q,W,U,V,SU,CV,XO,XEND)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,P,Q,S,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION CV(2,2),SU(2,2),U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,ALFA,ALFA2,AX,B,BETA,BETA2,BX,C,COSBX,COSCOS,
& COSHAX,CSHCSH,D,DETCU,DETCV,DETSU,DETSV,DETU,
& DETV,FAC,H,M,PRAD,RAD,RADA,RADB,RADM,SINBX,
& SINCOS,SINHAX,SINSIN,SNHCSH,SNHSNH,TR,TRUVA
C ..
C .. Local Arrays ..
DOUBLE PRECISION A(2,2),CUSVA(2,2),SUCVA(2,2),T(2,2),UUA(2,2),
& UVA(2,2),VUA(2,2),VVA(2,2)
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,EXP,SIGN,SIN,SQRT
C ..
C .. External Functions ..
DOUBLE PRECISION PHI
EXTERNAL PHI
C ..
D = 2.D0*SQRT(P* (Q-ELAM*W))
ALFA = 0.5D0*SQRT((S+D)/P)
BETA = 0.5D0*SQRT((-S+D)/P)
DETU = U(1,1)*U(2,2) - U(2,1)*U(1,2)
DETV = V(1,1)*V(2,2) - V(2,1)*V(1,2)
H = XEND - XO
AX = ABS(ALFA*H)
FAC = EXP(-AX)
IF (AX.GE.0.1D0) THEN
SINHAX = (1.D0-FAC*FAC)/2.D0
COSHAX = 1.D0 - SINHAX
SINHAX = SINHAX*SIGN(1.D0,H)/ALFA
ELSE
SINHAX = H*PHI(AX*AX)*FAC
COSHAX = (1.D0+FAC*FAC)/2.D0
END IF
BX = BETA*H
IF (ABS(BX).GE.0.1D0) THEN
SINBX = SIN(BX)*FAC/BETA
ELSE
SINBX = H*PHI(-BX*BX)*FAC
END IF
COSBX = COS(BX)*FAC
SNHCSH = SINHAX*COSHAX
SNHSNH = SINHAX*SINHAX
CSHCSH = COSHAX*COSHAX
SINCOS = SINBX*COSBX
SINSIN = SINBX*SINBX
COSCOS = COSBX*COSBX
ALFA2 = ALFA*ALFA
BETA2 = BETA*BETA
RAD = ALFA2 + BETA2
PRAD = P*RAD
RADM = ALFA2 - BETA2
RADA = 3.D0*ALFA2 - BETA2
RADB = 3.D0*BETA2 - ALFA2
C Compute SU.Adjoint(CV).exp(-AX):
SUCVA(1,1) = 0.5D0* (SNHCSH-SINCOS)/PRAD
SUCVA(1,2) = 0.25D0* (SNHSNH+SINSIN)/P
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = 0.5D0* (SNHCSH+SINCOS)/P
C Compute CU.Adjoint(SV).exp(-AX):
CUSVA(1,1) = 0.5D0* (RADA*SNHCSH-RADB*SINCOS)/P
CUSVA(1,2) = -0.25D0*RAD* (SNHSNH*BETA2-3.D0*CSHCSH+SINSIN*ALFA2+
& 3.D0*COSCOS)/P
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = 0.5D0*PRAD* (RADA*SNHCSH+RADB*SINCOS)
C Compute A = U.Adjoint(V)(xo):
A(1,1) = U(1,1)*V(2,2) - U(1,2)*V(2,1)
A(1,2) = U(1,2)*V(1,1) - U(1,1)*V(1,2)
A(2,1) = A(1,2)
A(2,2) = U(2,2)*V(1,1) - U(2,1)*V(1,2)
C Compute T = (CU.A.Adjoint(CV) + SU.Adjoint(A).Adjoint(SV))*exp(-AX):
T(1,1) = A(1,1)* (CSHCSH-SINSIN*BETA2) + A(1,2)* (SNHCSH+SINCOS) +
& A(2,2)* (CSHCSH-COSCOS)/RAD
T(1,2) = 0.5D0*RAD*A(1,1)* (SNHCSH-SINCOS) +
& 0.5D0*A(1,2)* (CSHCSH+COSCOS+BETA2*SNHSNH-ALFA2*SINSIN) +
& 0.5D0*A(2,2)* (SNHCSH+SINCOS)
T(2,1) = T(1,2)
T(2,2) = A(1,1)*RAD* (CSHCSH-COSCOS) +
& A(1,2)*RAD* (SNHCSH-SINCOS) +
& A(2,2)* (CSHCSH-SINSIN*BETA2)
C Compute U.Adjoint(V)(xend) = (CU.Uo+SU.Vo).Adjoint(SV.Uo+CV.Vo):
UVA(1,1) = DETU*CUSVA(1,1) + DETV*SUCVA(1,1) + T(1,1)
UVA(1,2) = DETU*CUSVA(1,2) + DETV*SUCVA(1,2) + T(1,2)
UVA(2,1) = DETU*CUSVA(2,1) + DETV*SUCVA(2,1) + T(2,1)
UVA(2,2) = DETU*CUSVA(2,2) + DETV*SUCVA(2,2) + T(2,2)
VUA(1,1) = UVA(2,2)
VUA(2,2) = UVA(1,1)
VUA(1,2) = -UVA(1,2)
VUA(2,1) = -UVA(2,1)
DETCU = 0.25D0*RADM* (SNHSNH+SINSIN) + 0.5D0* (CSHCSH+COSCOS)
DETSU = 0.25* (SNHSNH-SINSIN)/ (P*PRAD)
DETCV = DETCU
DETSV = 0.25D0*P*PRAD* (RADA*RADA*SNHSNH-RADB*RADB*SINSIN)
C Compute trace(CU.A.Adjoint(SU)).exp(-AX):
TR = 0.5D0* (A(1,1)* (SNHCSH+SINCOS)+A(1,2)* (SNHSNH+SINSIN))/P +
& 0.5D0*A(2,2)* (SNHCSH-SINCOS)/PRAD
UUA(1,1) = DETU*DETCU + DETV*DETSU + TR
UUA(1,2) = 0.D0
UUA(2,1) = 0.D0
UUA(2,2) = UUA(1,1)
C Compute trace(CV.Adjoint(A).Adjoint(SU)).exp(-AX):
TR = 0.5D0*P*RADA* (RAD* (A(1,1)*SNHCSH+A(1,2)*SNHSNH)+
& A(2,2)*SNHCSH) + 0.5D0*P*RADB*
& (RAD* (A(1,1)*SINCOS+A(1,2)*SINSIN)-A(2,2)*SINCOS)
VVA(1,1) = DETU*DETSV + DETV*DETCV + TR
VVA(1,2) = 0.D0
VVA(2,1) = 0.D0
VVA(2,2) = VVA(1,1)
TRUVA = UVA(1,1) + UVA(2,2)
AC = UUA(1,1) - VVA(1,1)
B = -TRUVA
C = UUA(1,1) + VVA(1,1)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
U(1,1) = UUA(1,1)*COS(M/2.D0) - UVA(1,1)*SIN(M/2.D0)
U(1,2) = -UVA(1,2)*SIN(M/2.D0)
U(2,1) = -UVA(2,1)*SIN(M/2.D0)
U(2,2) = UUA(2,2)*COS(M/2.D0) - UVA(2,2)*SIN(M/2.D0)
V(1,1) = VUA(1,1)*COS(M/2.D0) - VVA(1,1)*SIN(M/2.D0)
V(1,2) = VUA(1,2)*COS(M/2.D0)
V(2,1) = VUA(2,1)*COS(M/2.D0)
V(2,2) = VUA(2,2)*COS(M/2.D0) - VVA(2,2)*SIN(M/2.D0)
C That's done the very tricky renormalisation. Now do the SU and CV
C matrices which are the U and V matrices for Theta_{0}.
TRUVA = SUCVA(1,1) + SUCVA(2,2)
AC = (DETSU-DETCV)
B = -TRUVA
C = (DETSU+DETCV)
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
SU(1,1) = DETSU*COS(M/2.D0) + SUCVA(1,1)*SIN(M/2.D0)
SU(1,2) = -SUCVA(1,2)*SIN(M/2.D0)
SU(2,1) = -SUCVA(2,1)*SIN(M/2.D0)
SU(2,2) = DETSU*COS(M/2.D0) + SUCVA(2,2)*SIN(M/2.D0)
CV(1,1) = -SUCVA(2,2)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
CV(1,2) = -SUCVA(1,2)*COS(M/2.D0)
CV(2,1) = -SUCVA(2,1)*COS(M/2.D0)
CV(2,2) = -SUCVA(1,1)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE FUNSOL(XO,XEND,P,S,Q,W,ELAM,Y1,Y2,Y3,Y4,SCAL)
C .. Local Scalars ..
DOUBLE PRECISION ALFA,ALFA1,ALFAH,ALFASQ,BETA,BETA1,BETAH,BETASQ,
& COSBX,COSHAX,COSHHX,COSHNX,COSHX,COSNX,COSOMX,D,
& DISC,ETA,ETAH,ETASQ,FAC,FACETA,FACMU,FACNU,FACRH,
& H,MU,MUH,NU,NUH,NUSQ,OMEGA,OMEGA1,OMEGA2,OMH1,
& OMH2,OMSQ,RAD,RH,RHH,RHSQ,SINBX,SINHAX,SINHHX,
& SINHMX,SINHNX,SINHX,SINNX,SINOM1,SINOM2,SINOMX,
& SIZE,YTEMP2,YTEMP3
INTEGER I
C ..
C .. External Functions ..
C This subroutine finds the fundamental matrix for a fourth order
C Sturm-Liouville equation with constant coefficients
C
C py'''' - sy'' + qy = elam wy,
C
C on an interval [xo,xend]
C First measure the size of the coefficients.
DOUBLE PRECISION PHI
EXTERNAL PHI
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,P,Q,S,SCAL,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION Y1(4),Y2(4),Y3(4),Y4(4)
C ..
C .. External Subroutines ..
EXTERNAL FUNMAT
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,COS,EXP,SIN,SQRT
C ..
SCAL = 0.D0
H = XEND - XO
SIZE = ABS(H)* (ABS((ELAM*W-Q)/P)+ABS(S/P))
IF (ABS(SIZE).LT.0.02D0) THEN
C We can easily use the matrix exponential to get our fundamental
C solutions.
CALL FUNMAT(XO,XEND,P,S,Q,W,ELAM,Y1,Y2,Y3,Y4)
SCAL = 0.D0
ELSE
C We have to do the solutions the hard way. There are four cases
C to consider, each of which will have to be subdivided yet
C further.
ALFA1 = 0.5D0*S/P
DISC = S*S + 4.D0*P* (ELAM*W-Q)
BETA1 = 0.5D0*SQRT(ABS(DISC))/P
IF (DISC.LT.0.D0) THEN
C We have complex roots. We must compute their real and imaginary
C parts. This is tha case of large negative lambda.
D = 2.D0*SQRT(P* (Q-ELAM*W))
ALFA = 0.5D0*SQRT((S+D)/P)
BETA = 0.5D0*SQRT((D-S)/P)
ALFAH = ABS(ALFA*H)
SCAL = ALFAH
BETAH = BETA*H
FAC = EXP(-ALFAH)
COSHAX = (1.D0+FAC*FAC)/2.D0
IF (ALFAH.GT.0.1D0) THEN
SINHAX = (1.D0-FAC*FAC)/ (ALFAH+ALFAH)
ELSE
SINHAX = PHI(ALFAH*ALFAH)*FAC
END IF
SINHAX = H*SINHAX
COSBX = COS(BETAH)
IF (ABS(BETAH).GT.0.1D0) THEN
SINBX = SIN(BETAH)/BETAH
ELSE
SINBX = PHI(-BETAH*BETAH)
END IF
SINBX = H*SINBX
BETASQ = BETA*BETA
ALFASQ = ALFA*ALFA
Y1(1) = COSHAX*COSBX + (BETASQ-ALFASQ)*SINHAX*SINBX/2.D0
Y2(1) = (COSBX*SINHAX+COSHAX*SINBX)/2.D0
Y3(1) = (COSBX*SINHAX-COSHAX*SINBX)/
& (2.D0*P* (ALFASQ+BETASQ))
Y4(1) = SINHAX*SINBX/ (2.D0*P)
ELSE IF (DISC.GE.0.D0 .AND. ALFA1.GT.BETA1) THEN
C We have real roots only. This is the case of intermediate
C lambda and positive s.
NU = SQRT(ALFA1+BETA1)
MU = SQRT(ALFA1-BETA1)
ETA = 0.5D0* (NU+MU)
RH = 0.5D0* (NU-MU)
ETASQ = ETA*ETA
RHSQ = RH*RH
ETAH = ABS(ETA*H)
RHH = ABS(RH*H)
SCAL = ETAH + RHH
C This means that SCAL = ABS(NU*H)
FACETA = EXP(-ETAH)
FACRH = EXP(-RHH)
COSHNX = (1.D0+FACETA*FACETA)/2.D0
COSHHX = (1.D0+FACRH*FACRH)/2.D0
IF (ETAH.GT.0.1D0) THEN
SINHNX = (1.D0-FACETA*FACETA)/ (ETAH+ETAH)
ELSE
SINHNX = PHI(ETAH*ETAH)*FACETA
END IF
SINHNX = H*SINHNX
IF (RHH.GT.0.1D0) THEN
SINHHX = (1.D0-FACRH*FACRH)/ (RHH+RHH)
ELSE
SINHHX = PHI(RHH*RHH)*FACRH
END IF
SINHHX = H*SINHHX
Y1(1) = COSHNX*COSHHX - 0.5D0* (ETASQ+RHSQ)*SINHNX*SINHHX
Y2(1) = 0.5D0* (COSHNX*SINHHX+COSHHX*SINHNX)
Y4(1) = 0.5D0*SINHNX*SINHHX/P
IF (ETA*RH.GT.0.1D0) THEN
C This means that NU**2-MU**2 is not small
NUH = ABS(NU*H)
MUH = ABS(MU*H)
FACNU = EXP(-NUH)
FACMU = EXP(-MUH)
IF (NUH.GT.0.1D0) THEN
SINHNX = (1.D0-FACNU*FACNU)/ (NUH+NUH)
ELSE
SINHNX = PHI(NUH*NUH)*FACNU
END IF
SINHNX = H*SINHNX
IF (MUH.GT.0.1D0) THEN
SINHMX = (1.D0-FACMU*FACMU)/ (MUH+MUH)
ELSE
SINHMX = PHI(MUH*MUH)*FACMU
END IF
SINHMX = H*SINHMX*EXP(- (NUH-MUH))
Y3(1) = (SINHMX-SINHNX)/ (P* (NU**2-MU**2))
ELSE
Y3(1) = (COSHHX*SINHNX-COSHNX*SINHHX)/
& (2.D0*P* (ETASQ-RHSQ))
END IF
ELSE IF (DISC.GE.0.D0 .AND. ALFA1.LE.BETA1 .AND.
& ALFA1+BETA1.GE.0.D0) THEN
C We have two real roots and two complex roots. This is
C the case of large positive lambda.
NUSQ = ALFA1 + BETA1
OMSQ = BETA1 - ALFA1
NU = SQRT(NUSQ)
OMEGA = SQRT(OMSQ)
RAD = 2.D0*BETA1
NUH = ABS(NU*H)
SCAL = NUH
FAC = EXP(-NUH)
IF (NUH.GT.0.1D0) THEN
SINHNX = (1.D0-FAC*FAC)/ (NUH+NUH)
ELSE
SINHNX = PHI(NUH*NUH)*FAC
END IF
SINHNX = H*SINHNX
OMH1 = OMEGA*H
IF (ABS(OMH1).GT.0.1D0) THEN
SINOMX = SIN(OMH1)/OMH1
ELSE
SINOMX = PHI(-OMH1*OMH1)
END IF
SINOMX = H*SINOMX*FAC
COSOMX = COS(OMH1)*FAC
COSHNX = (1.D0+FAC*FAC)/2.D0
Y1(1) = (NUSQ*COSOMX+OMSQ*COSHNX)/RAD
Y2(1) = (NUSQ*SINHNX+OMSQ*SINOMX)/RAD
RAD = P*RAD
Y3(1) = (SINOMX-SINHNX)/RAD
Y4(1) = (COSHNX-COSOMX)/RAD
ELSE IF (DISC.GE.0.D0 .AND. ALFA1.LE.BETA1 .AND.
& ALFA1+BETA1.LE.0.D0) THEN
C This case is the case of purely complex roots; it is the
C case of intermediate lambda and negative s.
SCAL = 0.D0
OMEGA1 = SQRT(ABS(ALFA1+BETA1))
OMEGA2 = SQRT(ABS(ALFA1-BETA1))
ETA = 0.5D0* (OMEGA1+OMEGA2)
RH = 0.5D0* (OMEGA2-OMEGA1)
ETASQ = ETA*ETA
RHSQ = RH*RH
ETAH = ETA*H
RHH = RH*H
COSNX = COS(ETAH)
COSHX = COS(RHH)
IF (ABS(ETAH).GT.0.1D0) THEN
SINNX = SIN(ETAH)/ETAH
ELSE
SINNX = PHI(-ETAH*ETAH)
END IF
SINNX = H*SINNX
IF (ABS(RHH).GT.0.1D0) THEN
SINHX = SIN(RHH)/RHH
ELSE
SINHX = PHI(-RHH*RHH)
END IF
SINHX = H*SINHX
Y1(1) = COSNX*COSHX + (ETASQ+RHSQ)*SINNX*SINHX/2.D0
Y2(1) = (COSNX*SINHX+COSHX*SINNX)/2.D0
Y4(1) = SINNX*SINHX/ (2.D0*P)
IF (ETA*RH.GT.0.1D0) THEN
C OMEGA2**2-OMEGA1**2 is not small
OMH1 = OMEGA1*H
OMH2 = OMEGA2*H
IF (ABS(OMH1).GT.0.1D0) THEN
SINOM1 = SIN(OMH1)/OMH1
ELSE
SINOM1 = PHI(-OMH1*OMH1)
END IF
SINOM1 = H*SINOM1
IF (ABS(OMH2).GT.0.1D0) THEN
SINOM2 = SIN(OMH2)/OMH2
ELSE
SINOM2 = PHI(-OMH2*OMH2)
END IF
SINOM2 = H*SINOM2
Y3(1) = (SINOM2-SINOM1)/ (2.D0*P*BETA1)
ELSE
Y3(1) = - (COSHX*SINNX-COSNX*SINHX)/
& (2.D0*P* (ETASQ-RHSQ))
END IF
END IF
DO 10 I = 1,3
Y1(I+1) = (Q-ELAM*W)*Y3(I)
Y2(I+1) = Y1(I) + S*Y4(I)
Y3(I+1) = -Y4(I)
Y4(I+1) = Y2(I)/P
10 CONTINUE
YTEMP2 = Y1(3)
YTEMP3 = Y1(4)
Y1(4) = P*YTEMP2
Y1(3) = -P*YTEMP3 + S*Y1(2)
YTEMP2 = Y2(3)
YTEMP3 = Y2(4)
Y2(4) = P*YTEMP2
Y2(3) = -P*YTEMP3 + S*Y2(2)
YTEMP2 = Y3(3)
YTEMP3 = Y3(4)
Y3(4) = P*YTEMP2
Y3(3) = -P*YTEMP3 + S*Y3(2)
YTEMP2 = Y4(3)
YTEMP3 = Y4(4)
Y4(4) = P*YTEMP2
Y4(3) = -P*YTEMP3 + S*Y4(2)
END IF
RETURN
END
C -------------------------------------------------------------------
C This subroutine finds the exponential of a small matrix
C arising from a 4th order SL problem
SUBROUTINE FUNMAT(XO,XEND,P,S,Q,W,ELAM,Y1,Y2,Y3,Y4)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,P,Q,S,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION Y1(4),Y2(4),Y3(4),Y4(4)
C ..
C .. Local Scalars ..
DOUBLE PRECISION H,SCALAR
INTEGER I
C ..
C .. Local Arrays ..
DOUBLE PRECISION A(4,4),MNEST(4,4),STORE(4,4)
C ..
C .. External Subroutines ..
EXTERNAL F06QFF,F06QHF,NEST
C ..
C .. Intrinsic Functions ..
INTRINSIC REAL
C ..
H = XEND - XO
C Set up the matrix A.
C First set all elements to zero:
CALL F06QHF('G',4,4,0.D0,0.D0,A,4)
C Now fill in the non-zero elements:
A(1,2) = H
A(2,4) = H/P
A(3,1) = H* (Q-ELAM*W)
A(4,2) = H*S
A(4,3) = -H
C Now do the nested multiplication. Start by setting
C mnest to be the identity.
CALL F06QHF('G',4,4,0.D0,1.D0,MNEST,4)
DO 10 I = 9,1,-1
SCALAR = 1.D0/REAL(I)
CALL NEST(MNEST,4,4,SCALAR,A,4,STORE,4)
10 CONTINUE
C Now recover the matrix exponential into the vectors y1,y2,y3 and
C y4.
CALL F06QFF('G',4,1,MNEST(1,1),4,Y1,4)
CALL F06QFF('G',4,1,MNEST(1,2),4,Y2,4)
CALL F06QFF('G',4,1,MNEST(1,3),4,Y3,4)
CALL F06QFF('G',4,1,MNEST(1,4),4,Y4,4)
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE NEST(MNEST,IMNEST,N,SCALAR,FACTOR,IFACT,STORE,ISTOR)
C This routine performs the operation
C
C mnest = Identity + scalar.factor.mnest,
C
C where scalar is a real number, factor and mnest are nxn matrices,
C imnest is the first dimension of mnest as declared in the
C calling program and ifact is the first dimension of factor
C as declared in the calling program.
C
C .. Scalar Arguments ..
DOUBLE PRECISION SCALAR
INTEGER IFACT,IMNEST,ISTOR,N
C ..
C .. Array Arguments ..
DOUBLE PRECISION FACTOR(IFACT,N),MNEST(IMNEST,N),STORE(ISTOR,N)
C ..
C .. External Subroutines ..
EXTERNAL F06QFF,F06QHF,F06YAF
C ..
CALL F06QFF('G',N,N,MNEST,IMNEST,STORE,ISTOR)
CALL F06QHF('G',N,N,0.D0,1.D0,MNEST,IMNEST)
CALL F06YAF('N','N',N,N,N,SCALAR,FACTOR,IFACT,STORE,ISTOR,1.D0,
& MNEST,IMNEST)
RETURN
END
C ---------------------------------------------------------------------
SUBROUTINE SPTH4(U,V,OMEGA1,OMEGA2,ARG,LR,IFAIL)
C .. Scalar Arguments ..
DOUBLE PRECISION ARG,OMEGA1,OMEGA2
INTEGER IFAIL
CHARACTER LR
C ..
C .. Array Arguments ..
DOUBLE PRECISION U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION A,B,BLOC,C,DEFALT,DISC,PI,R,SGNA,SGNB,SGNC,TWOPI
C ..
C .. External Functions ..
DOUBLE PRECISION X01AAF
EXTERNAL X01AAF
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,MAX,MIN,SIGN,SQRT
C ..
C .. External Subroutines ..
C ..
IFAIL = 0
PI = X01AAF(0.D0)
TWOPI = 2.D0*PI
IF (LR.NE.'R') THEN
DEFALT = 0.D0
ELSE
DEFALT = TWOPI
END IF
A = V(1,1)*V(2,2) - V(1,2)*V(2,1)
B = U(1,2)*V(2,1) + U(2,1)*V(1,2) - U(1,1)*V(2,2) - U(2,2)*V(1,1)
C = U(1,1)*U(2,2) - U(1,2)*U(2,1)
R = MAX(ABS(A),ABS(B),ABS(C))
A = A/R
B = -B/R
BLOC = B
C = C/R
IF (B.EQ.0.D0) THEN
IF (C.EQ.0.D0) THEN
OMEGA1 = DEFALT
OMEGA2 = DEFALT
ELSE IF (A.EQ.0.D0) THEN
OMEGA1 = PI
OMEGA2 = OMEGA1
ELSE
OMEGA2 = ATAN2(-C,A)
IF (OMEGA1.LT.0.D0) OMEGA1 = OMEGA1 + TWOPI
IF (OMEGA2.LT.0.D0) OMEGA2 = OMEGA2 + TWOPI
END IF
ELSE
IF (A.EQ.0.D0) THEN
OMEGA1 = PI
IF (C.EQ.0.D0) THEN
OMEGA2 = DEFALT
ELSE
OMEGA2 = 2.D0*ATAN2(C,B)
IF (OMEGA2.LT.0.D0) OMEGA2 = OMEGA2 + TWOPI
END IF
ELSE
IF (C.EQ.0.D0) THEN
OMEGA1 = DEFALT
OMEGA2 = 2.D0*ATAN2(B,A)
IF (OMEGA2.LT.0.D0) OMEGA2 = OMEGA2 + TWOPI
ELSE
C We are now in the situation where A, B and C are all non-zero.
SGNA = SIGN(1.D0,A)
SGNB = SIGN(1.D0,B)
SGNC = SIGN(1.D0,C)
B = ABS(B)
DISC = SQRT(ABS(B*B-4.D0*A*C))
IF (SGNA*SGNC.GT.0.D0) THEN
OMEGA1 = 2.D0*ATAN2(SGNB*ABS(B-DISC),2.D0*A)
ELSE
OMEGA1 = 2.D0*ATAN2(-SGNB*ABS(B-DISC),2.D0*A)
END IF
OMEGA2 = 2.D0*ATAN2(SGNB*ABS(B+DISC),2.D0*A)
IF (OMEGA1.LT.0.D0) OMEGA1 = OMEGA1 + TWOPI
IF (OMEGA2.LT.0.D0) OMEGA2 = OMEGA2 + TWOPI
END IF
END IF
END IF
ARG = OMEGA1
OMEGA1 = MIN(OMEGA1,OMEGA2)
OMEGA2 = MAX(ARG,OMEGA2)
ARG = OMEGA1 + OMEGA2
B = BLOC
RETURN
END
C -----------------------------------------------------------------
SUBROUTINE SL4CNT(XEND,P,S,Q,W,ELAM,NZ)
C This routine counts the number of zeros of det(U(x)), where
C U is the 2x2 U-matrix of the Hamiltonian system associated
C with the 4th order SLDE (py'')''-(sy')'+qy = elam w y satisfying
C the initial condition U(0) = 0. The zeros counted are those
C which live in the range (0,xend] (or [xend,0) if xend<0).
C
C IMPORTANT NOTE: det(U(x)) is an EVEN function of x so we can
C always assume that xend>0.
C
C .. Scalar Arguments ..
C The values of xend, p, s, q, w, and elam are all supplied as
C input. The value of nz -- the number of zeros -- is output.
DOUBLE PRECISION ELAM,P,Q,S,W,XEND
INTEGER NZ
C ..
C .. Local Scalars ..
DOUBLE PRECISION CS,CSH,DET,DIFF,DISC,NU,OMEGA,OMEGA1,OMEGA2,PI,
& SDISC,SN,SNH,SUB,T,X
INTEGER ICASE,IZ
C ..
C .. External Functions ..
DOUBLE PRECISION PHI,X01AAF
EXTERNAL PHI,X01AAF
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,COS,EXP,INT,MOD,SIN,SQRT
C ..
PI = X01AAF(0.D0)
C The first thing to do is to classify the case that we are in.
DISC = S*S + 4.D0*P* (ELAM*W-Q)
SDISC = SQRT(ABS(DISC))
IF (DISC.LT.0.D0) THEN
C We are in the case of `large negative lambda'.
ICASE = 1
ELSE IF (DISC.GE.0.D0 .AND. S-SDISC.GE.0.D0) THEN
C This is the case of `Intermediate lambda and positive s'.
ICASE = 2
ELSE IF (DISC.GE.0.D0 .AND. S+SDISC.LE.0.D0) THEN
C This is the case of `Intermediate lambda and negative s'.
OMEGA1 = SQRT(0.5D0* (SDISC-S)/P)
OMEGA2 = SQRT(0.5D0* (-SDISC-S)/P)
ICASE = 3
ELSE IF (DISC.GE.0.D0 .AND. SDISC.GE.ABS(S)) THEN
C This is the case of `large positive lambda'.
NU = SQRT(0.5D0* (S+SDISC)/P)
OMEGA = SQRT(0.5D0* (SDISC-S)/P)
ICASE = 4
END IF
C Now we treat each case according to the theory in our paper.
C Cases 1 and 2 are a real scoosh, while 3 and 4 are much
C harder. Case 3 is the worst of the lot by a long way.
C Exploit the fact that det(U(x)) is an even function of x:
X = ABS(XEND)
GO TO (10,20,30,40) ICASE
C Case 1:
10 NZ = 0
RETURN
C Case 2:
20 NZ = 0
RETURN
C Case 3:
30 IF (OMEGA1.EQ.OMEGA2) THEN
NZ = 0
ELSE
IF (X*ABS(OMEGA1-OMEGA2).LE.2.D0*PI) THEN
NZ = 0
C A correction is necessary under certain circumstances
IF (X*ABS(OMEGA1+OMEGA2).GE.0.1D0) THEN
DIFF = ((OMEGA1-OMEGA2)**2)*
& (1.D0-COS((OMEGA1+OMEGA2)*X)) -
& ((OMEGA1+OMEGA2)**2)* (1.D0-
& COS((OMEGA1-OMEGA2)*X))
ELSE
DIFF = - (OMEGA1*OMEGA2* ((OMEGA1+OMEGA2)*X)**2*
& ((OMEGA1-OMEGA2)*X)**2)/6.D0
END IF
C IF (X*ABS(OMEGA1-OMEGA2).EQ.2.D0*PI .OR.
C & DIFF.LE.0.D0) NZ = 1
IF (X*ABS(OMEGA1-OMEGA2).EQ.2.D0*PI .OR.
& DIFF.GT.0.D0) NZ = 1
ELSE
IZ = INT(X*ABS(OMEGA1-OMEGA2)*0.5D0/PI)
NZ = 2*IZ - 1
C A correction is necessary under certain circumstances
IF (((OMEGA1-OMEGA2)**2)* (1.D0-COS((OMEGA1+
& OMEGA2)*X)).LE. ((OMEGA1+OMEGA2)**2)*
& (1.D0-COS((OMEGA1-OMEGA2)*X))) THEN
NZ = NZ + 1
ELSE IF (((OMEGA1-OMEGA2)**2)*
& (1.D0-COS((OMEGA1+OMEGA2)*X)).GE.
& ((OMEGA1+OMEGA2)**2)* (1.D0-COS((OMEGA1-
& OMEGA2)*X)) .AND. (OMEGA1+OMEGA2)*
& SIN((OMEGA1+OMEGA2)*X)* (1.D0-COS((OMEGA1-
& OMEGA2)*X)).GE. (OMEGA1-OMEGA2)*
& SIN((OMEGA1-OMEGA2)*X)* (1.D0-COS((OMEGA1+
& OMEGA2)*X))) THEN
NZ = NZ + 2
END IF
END IF
END IF
RETURN
C Case 4:
40 IF (OMEGA.EQ.0.D0 .OR. X*OMEGA.LE.PI) THEN
C There are no zeros in this case
NZ = 0
C ???
T = NU*X
IF (T.GT.0.1D0) THEN
SNH = 0.5D0* (1.D0-EXP(-2.D0*T))/T
CSH = 0.5D0* (1.D0+EXP(-2.D0*T))
SUB = EXP(-T)
ELSE
SNH = PHI(T*T)
CSH = 0.5D0* (EXP(T)+EXP(-T))
SUB = 1.D0
END IF
T = OMEGA*X
CS = COS(T)
IF (T.GT.0.1D0) THEN
SN = SIN(T)/T
ELSE
SN = PHI(-T*T)
END IF
DET = CS*CSH - SUB - 0.5D0* (NU**2-OMEGA**2)*SN*SNH
C ???
ELSE
NZ = INT(X*OMEGA/PI) - 1
C We may have to increment nz by 1 to get the zero-count right
T = NU*X
IF (T.GT.0.1D0) THEN
SNH = 0.5D0* (1.D0-EXP(-2.D0*T))/T
CSH = 0.5D0* (1.D0+EXP(-2.D0*T))
SUB = EXP(-T)
ELSE
SNH = PHI(T*T)
CSH = 0.5D0* (EXP(T)+EXP(-T))
SUB = 1.D0
END IF
T = OMEGA*X
CS = COS(T)
IF (T.GT.0.1D0) THEN
SN = SIN(T)/T
ELSE
SN = PHI(-T*T)
END IF
DET = CS*CSH - SUB - 0.5D0* (NU**2-OMEGA**2)*SN*SNH
C If the determinant has not changed sign and the number of zeros in
C previous intervals is odd then there must be a further zero to make
C the total number of zeros even.
IF (DET.LE.0.D0 .AND. MOD(NZ,2).NE.0) THEN
NZ = NZ + 1
ELSE IF (DET.GE.0.D0 .AND. MOD(NZ,2).EQ.0) THEN
NZ = NZ + 1
END IF
C If the determinant has changed sign and the number of zeros in
C previous intervals is even then there must be a further zero to make
C the total number of zeros odd.
END IF
RETURN
END
C =====================================================================
C ============== SUBROUTINES FOR EIGENFUNCTION APPROXIMATION ==========
C =====================================================================
SUBROUTINE EFN(XMESH,P,S,Q,W,ELAM,X,YSTOR1,IMATCH,NMESH,Y,HOT)
C
C This subroutines takes as input the arrays containing the nodal values
C of the eigenfunction, and uses Hermite interpolation to determine the
C values of the eigenfunction and of its quasiderivatives at internodal
C points.
C
C The parameter HOT should be set to .FALSE. on the first call. It will
C be reset to .TRUE. by EFN; on subsequent calls for computation of the
C same eigenfunction, HOT should not be changed by the user. This allows
C the routine to reduce the amount of searching required to find a mesh
C interval containing the point at which the eigenfunction is to be
C evaluated.
C
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,X
INTEGER IMATCH,NMESH
LOGICAL HOT
C ..
C .. Array Arguments ..
DOUBLE PRECISION P(1:NMESH),Q(1:NMESH),S(1:NMESH),W(1:NMESH),
& XMESH(0:NMESH),Y(1:4),YSTOR1(4,0:NMESH)
C ..
C .. Local Scalars ..
DOUBLE PRECISION EPS,LSTVLU,X1,X2,XLOC
INTEGER I,J,K,L,M
C ..
C .. External Subroutines ..
EXTERNAL AEFP
C ..
C .. External Functions ..
DOUBLE PRECISION X02AJF
EXTERNAL X02AJF
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,MAX,MIN
C ..
C .. Save statement ..
SAVE EPS,LSTVLU,M
C ..
IF (.NOT.HOT) THEN
J = 1
K = NMESH
L = 1
LSTVLU = XMESH(0)
EPS = 100.D0*X02AJF(0.D0)
HOT = .TRUE.
GO TO 10
END IF
IF (X.GT.LSTVLU) THEN
J = M
K = NMESH
L = 1
ELSE
J = M
K = 1
L = -1
END IF
10 XLOC = MIN(X,XMESH(NMESH))
XLOC = MAX(XLOC,XMESH(0))
DO 20 I = J,K,L
X1 = XMESH(I-1)
X2 = XMESH(I)
IF (XLOC.LT.X2 .AND. XLOC.GE.X1) THEN
IF (ABS(XLOC-X2).LT.EPS) THEN
C Just assign the value at X2
Y(4) = YSTOR1(4,I)
Y(3) = YSTOR1(3,I)
Y(2) = YSTOR1(2,I)
Y(1) = YSTOR1(1,I)
ELSE IF (ABS(XLOC-X1).LT.EPS) THEN
C Just assign the value at X1
Y(4) = YSTOR1(4,I-1)
Y(3) = YSTOR1(3,I-1)
Y(2) = YSTOR1(2,I-1)
Y(1) = YSTOR1(1,I-1)
ELSE
C Do Hermite interpolation
CALL AEFP(X1,X2,XLOC,YSTOR1(4,I-1),
& (S(I)*YSTOR1(2,I-1)-YSTOR1(3,I-1)),
& YSTOR1(4,I), (S(I)*YSTOR1(2,I)-YSTOR1(3,I)),
& Y(4))
CALL AEFP(X1,X2,XLOC,YSTOR1(1,I-1),YSTOR1(2,I-1),
& YSTOR1(1,I),YSTOR1(2,I),Y(1))
CALL AEFP(X1,X2,XLOC,YSTOR1(3,I-1),
& (Q(I)-ELAM*W(I))*YSTOR1(1,I-1),YSTOR1(3,I),
& (Q(I)-ELAM*W(I))*YSTOR1(1,I-1),Y(3))
CALL AEFP(X1,X2,XLOC,YSTOR1(2,I-1),
& (YSTOR1(4,I-1)/P(I)),YSTOR1(2,I),
& (YSTOR1(4,I)/P(I)),Y(2))
END IF
M = I
LSTVLU = XLOC
GO TO 30
END IF
20 CONTINUE
30 CONTINUE
RETURN
END
C ---------------------------------------------------------------------
SUBROUTINE AEFP(XO,XEND,INTP,FXO,FDXO,FXEND,FDXEND,Y)
C .. Scalar Arguments ..
DOUBLE PRECISION FDXEND,FDXO,FXEND,FXO,INTP,XEND,XO,Y
C ..
C .. Local Scalars ..
DOUBLE PRECISION H,H2,H3,XA,XA2,XB,XB2
C ..
H = XEND - XO
H2 = H*H
H3 = H*H2
XA = INTP - XEND
XB = INTP - XO
XA2 = XA*XA
XB2 = XB*XB
Y = FXO* ((XA2/H2)+ ((2.D0*XB*XA2)/H3)) +
& FXEND* ((XB2/H2)- ((2.D0*XA*XB2)/H3)) + FDXO* (XB*XA2/H2) +
& FDXEND* (XA*XB2/H2)
RETURN
END
C ----------------------------------------------------------------------
SUBROUTINE SL4EFN(X,Y1,Y2,HOT,WORK,IWORK,NMESH,IMATCH,ELAM,MULTI,
& WORKS,YSTOR1,YSTOR2,SL4BCS,IFAIL)
C This subroutine computes eigenfunctions of fourth-order
C Sturm-Liouville problems.
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,X
INTEGER IFAIL,IMATCH,IWORK,MULTI,NMESH
LOGICAL HOT
C ..
C .. Array Arguments ..
DOUBLE PRECISION WORK(0:IWORK,1:6),WORKS(19*NMESH+48),Y1(1:4),
& Y2(1:4),YSTOR1(1:4,0:NMESH),YSTOR2(1:4,0:NMESH)
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4BCS
C ..
C .. Local Scalars ..
DOUBLE PRECISION ARGL,EPS,FAC,NUX,NUX2
INTEGER FMTCH1,FNMSH1,I,IAL,IAR,IDL1,IDL2,IDR1,IDR2,IFO,IMTCH1,
& INDEX,INDEX2,INL,INR,IRL,IRR,IUL,IUR,IV,IVL,IVR,IW,IX,IZL,
& IZR,IZVL,IZVR,J,NMSH1,NREC,TMTCH1,TNMSH1
LOGICAL PRN
CHARACTER*6 SRNAME
C ..
C .. Local Arrays ..
DOUBLE PRECISION RINV(2,2),U(2,2),ULOC(2,2),V(2,2),VLOC(2,2),
& YLOC(4)
CHARACTER*80 REC(2)
C ..
C .. External Functions ..
DOUBLE PRECISION X02AJF
INTEGER P01ABF
EXTERNAL X02AJF,P01ABF
C ..
C .. External Subroutines ..
EXTERNAL F06PAF,GETVEC,THTNT1
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,EXP,MAX,MIN
C ..
C .. Save statement ..
SAVE IMTCH1,TMTCH1,FMTCH1,NMSH1,TNMSH1,FNMSH1,INL,INR,IRL,IRR,IUL,
& IUR,IVL,IVR,IV,IW,IZL,IZR,IZVL,IZVR,IAL,IAR,IDL1,IDL2,IDR1,
& IDR2
C ..
SRNAME = 'SL4EFN'
IFO = IFAIL
IF (ABS(IFO).GT.1) IFO = 0
IFAIL = 0
EPS = 1000.D0*X02AJF(0.D0)
IF (.NOT.HOT) THEN
IMTCH1 = IMATCH + 1
TMTCH1 = 2*IMTCH1
FMTCH1 = 4*IMTCH1
NMSH1 = NMESH - IMATCH + 1
TNMSH1 = 2*NMSH1
FNMSH1 = 4*NMSH1
INL = 1
INR = INL + IMTCH1
IRL = INR + NMSH1
IRR = IRL + FMTCH1
IUL = IRR + FNMSH1
IUR = IUL + FMTCH1
IVL = IUR + FNMSH1
IVR = IVL + FMTCH1
IV = IVR + FNMSH1
IW = IV + TNMSH1
IZL = IW + TMTCH1
IZR = IZL + TMTCH1
IZVL = IZR + TNMSH1
IZVR = IZVL + TMTCH1
IAL = IZVR + TNMSH1
IAR = IAL + 1
IDL1 = IAR + 1
IDL2 = IDL1 + 2
IDR1 = IDL2 + 2
IDR2 = IDR1 + 2
C We must start by evaluating the eigenfunction at all the mesh-points.
C NOTE: the unextrapolated value of the eigenvalue is hidden in
C WORK(0,3).
C
CALL GETVEC(WORK(0,1),NMESH,IMATCH,WORK(1,2),WORK(1,3),
& WORK(1,4),WORK(1,5),WORK(0,3),MULTI,WORKS,YSTOR1,
& YSTOR2,SL4BCS,IFAIL)
IF (IFAIL.NE.0) GO TO 250
C We have VSTORE and WSTORE in WORKS(IV..) and WORKS(IW..). If the
C eigenfunction in question has multiplicity greater than 1, then we
C must recover the second copy of VSTOR and WSTOR.
IF (MULTI.EQ.2) THEN
WORKS(IZVL+2*IMATCH) = WORKS(IDL1)
WORKS(IZVL+2*IMATCH+1) = WORKS(IDL1+1)
DO 10 I = IMATCH,1,-1
CALL F06PAF('N',2,2,1.D0,WORKS(IRL+4*I),2,
& WORKS(IZVL+2*I),1,0.D0,
& WORKS(IZVL+2* (I-1)),1)
10 CONTINUE
WORKS(IZVR) = WORKS(IDR1)
WORKS(IZVR+1) = WORKS(IDR1+1)
DO 20 I = IMATCH,NMESH
CALL F06PAF('N',2,2,1.D0,WORKS(IRR+4* (I-IMATCH)),2,
& WORKS(IZVR+2* (I-IMATCH)),1,0.D0,
& WORKS(IZVR+2* (I-IMATCH+1)),1)
20 CONTINUE
END IF
END IF
C We are now ready to get the eigenfunction value at the point X
C specified by the user.
IF (X.LT.WORK(0,1)-EPS .OR. X.GT.WORK(NMESH,1)+EPS) THEN
IFAIL = 1
GO TO 250
END IF
X = MAX(X,WORK(0,1))
X = MIN(X,WORK(NMESH,1))
IX = 0
30 IF (WORK(IX+1,1).GE.X .OR. IX+1.GE.NMESH) GO TO 40
IX = IX + 1
GO TO 30
C Now we have found an interval [WORK(IX,1),WORK(IX+1,1)] containing
C the point X. We integrate over the interval to find the value of the
C eigenfunction at X from the endpoint values.
40 IF (X.LT.WORK(IX,1)+1.D2*X02AJF(0.D0)) THEN
DO 50 I = 1,4
Y1(I) = YSTOR1(I,IX)
50 CONTINUE
IF (MULTI.EQ.2) THEN
DO 60 I = 1,4
Y2(I) = YSTOR2(I,IX)
60 CONTINUE
END IF
ELSE IF (X.GT.WORK(IX+1,1)-1.D2*X02AJF(0.D0)) THEN
DO 70 I = 1,4
Y1(I) = YSTOR1(I,IX+1)
70 CONTINUE
IF (MULTI.EQ.2) THEN
DO 80 I = 1,4
Y2(I) = YSTOR2(I,IX+1)
80 CONTINUE
END IF
ELSE
C In this case we are not close to one of the endpoints so we must do
C some work. First, we need the integers which indicate the position
C in the storage array WORKS of various essential animals:
PRN = .FALSE.
ARGL = 0.D0
NUX = 0.D0
NUX2 = 0.D0
IFAIL = 0
IF (IX.LT.IMATCH) THEN
C We get our eigenfunction value with a forward integration.
INDEX = IUL + 4*IX
INDEX2 = IVL + 4*IX
DO 100 J = 1,2
DO 90 I = 1,2
U(I,J) = WORKS(INDEX)
V(I,J) = WORKS(INDEX2)
INDEX = INDEX + 1
INDEX2 = INDEX2 + 1
90 CONTINUE
100 CONTINUE
CALL THTNT1(U,V,RINV,ARGL,NUX,WORK(IX+1,2),WORK(IX+1,3),
& WORK(IX+1,4),WORK(IX+1,5),WORK(0,3),
& WORK(IX,1),X,PRN,IFAIL)
IF (IFAIL.NE.0) GO TO 250
DO 120 J = 1,2
DO 110 I = 1,2
ULOC(I,J) = U(I,J)
VLOC(I,J) = V(I,J)
110 CONTINUE
120 CONTINUE
CALL THTNT1(ULOC,VLOC,RINV,ARGL,NUX2,WORK(IX+1,2),
& WORK(IX+1,3),WORK(IX+1,4),WORK(IX+1,5),
& WORK(0,3),X,WORK(IX+1,1),PRN,IFAIL)
IF (IFAIL.NE.0) GO TO 250
C
C We now have all the information that we need in order to recover
C the eigenfunction. We use the formula
C
C / y(x) \ = U.RINV.WSTOR_{IX+1}.
C \ y'(x) / ALFAL.EXP(-NUX2-SUM_{I=IX+2}^{IMATCH}NLSTOR(I))
C
C / -(py''(x))'+sy'(x) \ = V.RINV.WSTOR_{IX+1}.ALFAL.
C \ py''(x) / EXP(-NUX2-SUM_{I=IX+2}^{IMATCH}NLSTOR(I))
C
C valid for x < xmesh(IMATCH). WSTOR is in WORKS(IW+2*(IX+1)). ALFAL is
C in WORKS(IAL). NLSTOR(I) is in WORKS(INL+I).
C
IF (MULTI.NE.2) THEN
CALL F06PAF('N',2,2,1.D0,RINV,2,WORKS(IW+2* (IX+1)),1,
& 0.D0,YLOC(1),1)
CALL F06PAF('N',2,2,1.D0,U,2,YLOC(1),1,0.D0,Y1(1),1)
CALL F06PAF('N',2,2,1.D0,V,2,YLOC(1),1,0.D0,Y1(3),1)
C Now do the scaling.
NUX = NUX2
IF (IX+1.LT.IMATCH) THEN
DO 130 I = IX + 2,IMATCH
NUX = NUX + WORKS(INL+I)
130 CONTINUE
END IF
FAC = EXP(-NUX)*WORKS(IAL)
DO 140 I = 1,4
Y1(I) = Y1(I)*FAC
140 CONTINUE
ELSE
C We must also get the second eigenfunction.
CALL F06PAF('N',2,2,1.D0,RINV,2,WORKS(IW+2* (IX+1)),1,
& 0.D0,YLOC(1),1)
CALL F06PAF('N',2,2,1.D0,U,2,YLOC(1),1,0.D0,Y2(1),1)
CALL F06PAF('N',2,2,1.D0,V,2,YLOC(1),1,0.D0,Y2(3),1)
CALL F06PAF('N',2,2,1.D0,RINV,2,WORKS(IZVL+2* (IX+1)),
& 1,0.D0,YLOC(1),1)
CALL F06PAF('N',2,2,1.D0,U,2,YLOC(1),1,0.D0,Y1(1),1)
CALL F06PAF('N',2,2,1.D0,V,2,YLOC(1),1,0.D0,Y1(3),1)
C Now do the scaling.
NUX = NUX2
IF (IX+1.LT.IMATCH) THEN
DO 150 I = IX + 2,IMATCH
NUX = NUX + WORKS(INL+I)
150 CONTINUE
END IF
FAC = EXP(-NUX)*WORKS(IAL)
DO 160 I = 1,4
Y1(I) = Y1(I)*FAC
Y2(I) = Y2(I)*FAC
160 CONTINUE
END IF
ELSE
C We get our eigenfunction value with a backward integration.
INDEX = IUR + 4* (IX+1-IMATCH)
INDEX2 = IVR + 4* (IX+1-IMATCH)
DO 180 J = 1,2
DO 170 I = 1,2
U(I,J) = WORKS(INDEX)
V(I,J) = WORKS(INDEX2)
INDEX = INDEX + 1
INDEX2 = INDEX2 + 1
170 CONTINUE
180 CONTINUE
CALL THTNT1(U,V,RINV,ARGL,NUX,WORK(IX+1,2),WORK(IX+1,3),
& WORK(IX+1,4),WORK(IX+1,5),WORK(0,3),
& WORK(IX+1,1),X,PRN,IFAIL)
IF (IFAIL.NE.0) GO TO 250
DO 200 J = 1,2
DO 190 I = 1,2
ULOC(I,J) = U(I,J)
VLOC(I,J) = V(I,J)
190 CONTINUE
200 CONTINUE
CALL THTNT1(ULOC,VLOC,RINV,ARGL,NUX2,WORK(IX+1,2),
& WORK(IX+1,3),WORK(IX+1,4),WORK(IX+1,5),
& WORK(0,3),X,WORK(IX,1),PRN,IFAIL)
IF (IFAIL.NE.0) GO TO 250
C
C We now have all the information that we need in order to recover
C the eigenfunction. We use the formula
C
C / y(x) \ = U.RINV.VSTOR_{IX}.ALFAR.
C \ y'(x) / EXP(-NUX2-SUM_{I=IMATCH}^{IX-1}NRSTOR(I))
C
C / -(py''(x))'+sy'(x) \ = V.RINV.VSTOR_{IX}.ALFAR.
C \ py''(x) / EXP(-NUX2-SUM_{I=IMATCH}^{IX-1}NRSTOR(I))
C
C valid for x > xmesh(IMATCH). VSTOR is in WORKS(IV+2*(IX-IMATCH)).
C ALFAR is in WORKS(IAR). NRSTOR(I) is in WORKS(INR+I-1).
C
IF (MULTI.NE.2) THEN
CALL F06PAF('N',2,2,1.D0,RINV,2,
& WORKS(IV+2* (IX-IMATCH)),1,0.D0,YLOC(1),1)
CALL F06PAF('N',2,2,1.D0,U,2,YLOC(1),1,0.D0,Y1(1),1)
CALL F06PAF('N',2,2,1.D0,V,2,YLOC(1),1,0.D0,Y1(3),1)
C Now do the scaling.
NUX = NUX2
IF (IX.GT.IMATCH) THEN
DO 210 I = 0,IX - IMATCH - 1
NUX = NUX + WORKS(INR+I)
210 CONTINUE
END IF
FAC = EXP(-NUX)*WORKS(IAR)
DO 220 I = 1,4
Y1(I) = Y1(I)*FAC
220 CONTINUE
ELSE
C We must also get the second eigenfunction.
CALL F06PAF('N',2,2,1.D0,RINV,2,
& WORKS(IV+2* (IX-IMATCH)),1,0.D0,YLOC(1),1)
CALL F06PAF('N',2,2,1.D0,U,2,YLOC(1),1,0.D0,Y2(1),1)
CALL F06PAF('N',2,2,1.D0,V,2,YLOC(1),1,0.D0,Y2(3),1)
CALL F06PAF('N',2,2,1.D0,RINV,2,
& WORKS(IZVR+2* (IX-IMATCH)),1,0.D0,YLOC(1),
& 1)
CALL F06PAF('N',2,2,1.D0,U,2,YLOC(1),1,0.D0,Y1(1),1)
CALL F06PAF('N',2,2,1.D0,V,2,YLOC(1),1,0.D0,Y1(3),1)
C Now do the scaling.
NUX = NUX2
IF (IX.GT.IMATCH) THEN
DO 230 I = 0,IX - IMATCH - 1
NUX = NUX + WORKS(INR+I)
230 CONTINUE
END IF
FAC = EXP(-NUX)*WORKS(IAR)
DO 240 I = 1,4
Y1(I) = Y1(I)*FAC
Y2(I) = Y2(I)*FAC
240 CONTINUE
END IF
END IF
END IF
HOT = .TRUE.
250 IF (IFAIL.EQ.1) THEN
WRITE (REC,FMT=9000)
9000 FORMAT (' ** The X which you specified was out of range')
NREC = 1
ELSE IF (ABS(IFAIL).GT.1) THEN
IFAIL = 2
WRITE (REC,FMT=9010)
9010 FORMAT (' ** Some Theta matrices have arisen whose',/,
& ' ** eigenvalues cannot be accurately determined')
NREC = 2
END IF
IFAIL = P01ABF(IFO,IFAIL,SRNAME,NREC,REC)
RETURN
END
C ----------------------------------------------------------------------
SUBROUTINE GETVCN(WORK,IWORK,NMESH,IMATCH,ELAM,MULTI,WORKS,YSTOR1,
& YSTOR2,SL4BCS,IFAIL)
C INPUT VARIABLES
C ---------------
C WORK is input from a previous call to SLEUTH, and contains the
C information about the mesh and the coefficient functions PP, SP, QP
C and WP. It must be of length (0:IWORK,1:6), where IWORK is the same
C as on the call to SLEUTH.
C
C NMESH is the value returned by SLEUTH, as is IMATCH.
C ELAM is the eigenvalue approximation returned by SLEUTH
C
C OUTPUT VARIABLES
C ----------------
C MULTI is the multiplicity of ELAM as an eigenvalue of the problem;
C MULTI = 1 or 2.
C
C YSTOR1 is aan array of length (1:4,0:NMESH) such that
C YSTOR1(1,i) contains y(x(i)) for i = 0,..,NMESH
C YSTOR1(2,i) contains y'(x(i)) for i = 0,..,NMESH
C YSTOR1(3,i) contains -(py'')'(x(i))+sy'(x(i)) for i = 0,..,NMESH
C YSTOR1(4,i) contains py''(x(i)) for i = 0,..,NMESH
C
C YSTOR2 contains the same information for the second eigenfunction
C if MULTI = 2.
C
C WORKSPACE
C ---------
C WORKS is an array of length (1:19*NMESH+48). It is used as workspace.
C
C ROUTINES SUPPLIED BY THE USER
C -----------------------------
C SL4BCS is the same routine as was used by SLEUTH.
C
C
C Remark: WORK(0,3) contains, hidden, the last un-extrapolated
C eigenvalue approximation to be calculated.
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM
INTEGER IFAIL,IMATCH,IWORK,MULTI,NMESH
C ..
C .. Array Arguments ..
DOUBLE PRECISION WORK(0:IWORK,1:6),WORKS(19*NMESH+48),
& YSTOR1(1:4,0:NMESH),YSTOR2(1:4,0:NMESH)
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4BCS
C ..
C .. External Subroutines ..
EXTERNAL GETVEC
C ..
CALL GETVEC(WORK(0,1),NMESH,IMATCH,WORK(1,2),WORK(1,3),WORK(1,4),
& WORK(1,5),WORK(0,3),MULTI,WORKS,YSTOR1,YSTOR2,SL4BCS,
& IFAIL)
RETURN
END
C ---------------------------------------------------------------------
SUBROUTINE GETVEC(XMESH,NMESH,IMATCH,PP,SP,QP,WP,ELAM,MULTI,WORKS,
& YSTOR1,YSTOR2,SL4BCS,IFAIL)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM
INTEGER IFAIL,IMATCH,MULTI,NMESH
C ..
C .. Array Arguments ..
C WORKS must be of length at least 19*NMESH+48
DOUBLE PRECISION PP(1:NMESH),QP(1:NMESH),SP(1:NMESH),
& WORKS(19*NMESH+48),WP(1:NMESH),XMESH(0:NMESH),
& YSTOR1(4,0:NMESH),YSTOR2(4,0:NMESH)
C ..
C .. External Subroutines ..
EXTERNAL GETVC1
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4BCS
C ..
C IMPORTANT NOTE: The array WORKS must be of length at least
C 19*(NMESH+2)+10
C .. Local Scalars ..
INTEGER FMTCH1,FNMSH1,IAL,IAR,IDL1,IDL2,IDR1,IDR2,IMTCH1,INL,INR,
& IRL,IRR,IUL,IUR,IV,IVL,IVR,IW,IZL,IZR,IZVL,IZVR,NMSH1,
& TMTCH1,TNMSH1
C ..
IFAIL = 0
IMTCH1 = IMATCH + 1
TMTCH1 = 2*IMTCH1
FMTCH1 = 4*IMTCH1
NMSH1 = NMESH - IMATCH + 1
TNMSH1 = 2*NMSH1
FNMSH1 = 4*NMSH1
INL = 1
INR = INL + IMTCH1
IRL = INR + NMSH1
IRR = IRL + FMTCH1
IUL = IRR + FNMSH1
IUR = IUL + FMTCH1
IVL = IUR + FNMSH1
IVR = IVL + FMTCH1
IV = IVR + FNMSH1
IW = IV + TNMSH1
IZL = IW + TMTCH1
IZR = IZL + TMTCH1
IZVL = IZR + TNMSH1
IZVR = IZVL + TMTCH1
IAL = IZVR + TNMSH1
IAR = IAL + 1
IDL1 = IAR + 1
IDL2 = IDL1 + 2
IDR1 = IDL2 + 2
IDR2 = IDR1 + 2
CALL GETVC1(XMESH,NMESH,IMATCH,PP,SP,QP,WP,ELAM,MULTI,WORKS(INL),
& WORKS(INR),WORKS(IRL),WORKS(IRR),WORKS(IUL),
& WORKS(IUR),WORKS(IVL),WORKS(IVR),WORKS(IV),WORKS(IW),
& WORKS(IZL),WORKS(IZR),WORKS(IZVL),WORKS(IZVR),YSTOR1,
& YSTOR2,WORKS(IAL),WORKS(IAR),WORKS(IDL1),WORKS(IDL2),
& WORKS(IDR1),WORKS(IDR2),SL4BCS,IFAIL)
RETURN
END
C ---------------------------------------------------------------------
SUBROUTINE GETVC1(XMESH,NMESH,IMATCH,PP,SP,QP,WP,ELAM,MULTI,
& NLSTOR,NRSTOR,RLSTOR,RRSTOR,ULSTOR,URSTOR,
& VLSTOR,VRSTOR,VSTORE,WSTORE,ZLSTOR,ZRSTOR,
& ZVLSTR,ZVRSTR,YSTOR1,YSTOR2,ALFAL,ALFAR,DL1,DL2,
& DR1,DR2,SL4BCS,IFAIL)
C .. Scalar Arguments ..
DOUBLE PRECISION ALFAL,ALFAR,ELAM
INTEGER IFAIL,IMATCH,MULTI,NMESH
C ..
C .. Array Arguments ..
DOUBLE PRECISION DL1(2),DL2(2),DR1(2),DR2(2),NLSTOR(0:IMATCH),
& NRSTOR(IMATCH:NMESH),PP(1:NMESH),QP(1:NMESH),
& RLSTOR(2,2,0:IMATCH),RRSTOR(2,2,IMATCH:NMESH),
& SP(1:NMESH),ULSTOR(2,2,0:IMATCH),
& URSTOR(2,2,IMATCH:NMESH),VLSTOR(2,2,0:IMATCH),
& VRSTOR(2,2,IMATCH:NMESH),VSTORE(2,IMATCH:NMESH),
& WP(1:NMESH),WSTORE(2,0:IMATCH),XMESH(0:NMESH),
& YSTOR1(4,0:NMESH),YSTOR2(4,0:NMESH),
& ZLSTOR(2,0:IMATCH),ZRSTOR(2,IMATCH:NMESH),
& ZVLSTR(2,0:IMATCH),ZVRSTR(2,IMATCH:NMESH)
C ..
C .. Local Scalars ..
DOUBLE PRECISION ARGL,ARGR,BETAL,BETAR,DUMMY2,EXCONT,MAXUC,NUX,
& RATIO,XEND,XO
INTEGER I,J,JJ,K
LOGICAL ISING,LAST,PRN
C ..
C .. Local Arrays ..
DOUBLE PRECISION DUMMY(4),RINV(2,2),UC(2,2),UL(2,2),UR(2,2),
& VL(2,2),VR(2,2)
C ..
C .. External Functions ..
DOUBLE PRECISION X02AJF
EXTERNAL X02AJF
C ..
C .. External Subroutines ..
EXTERNAL F06PAF,F06QFF,F06QHF,F06YAF,THTNT1
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4BCS
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,EXP,LOG,MAX,SQRT
C ..
IFAIL = 0
PRN = .FALSE.
C Get the boundary condition at the left-hand endpoint:
CALL SL4BCS(0,ISING,XMESH(0),UL,VL,ELAM)
CALL F06QFF('G',2,2,UL,2,ULSTOR(1,1,0),2)
CALL F06QFF('G',2,2,VL,2,VLSTOR(1,1,0),2)
CALL F06QHF('G',2,2,0.D0,1.D0,RLSTOR(1,1,0),2)
NLSTOR(0) = 0.D0
C Now step across the mesh from XMESH(0) to XMESH(IMATCH),
C carrying the matrices UL, VL such that Theta = (VL+iUL).inv(VL-iUL)
C ARGL is not important so we set it (arbitarily) to zero.
ARGL = 0.D0
DO 30 I = 1,IMATCH
XO = XMESH(I-1)
XEND = XMESH(I)
C The routine THTNT1 does the stepping across the interval
C [XMESH(I-1),XMESH(I)] on which the coefficients P, S, Q and
C W have the constant values PP(I), SP(I), QP(I) and WP(I).
C The eigenparameter has the value ELAM. PRN is a logical
C variable which controls printing of intermediate results
C and IFAIL is an integer which is used to flag errors.
CALL THTNT1(UL,VL,RINV,ARGL,NUX,PP(I),SP(I),QP(I),WP(I),ELAM,
& XO,XEND,PRN,IFAIL)
IF (IFAIL.NE.0) RETURN
DO 20 K = 1,2
DO 10 J = 1,2
ULSTOR(J,K,I) = UL(J,K)
VLSTOR(J,K,I) = VL(J,K)
RLSTOR(J,K,I) = RINV(J,K)
10 CONTINUE
20 CONTINUE
NLSTOR(I) = NUX
30 CONTINUE
C Get the boundary condition at the right-hand endpoint:
CALL SL4BCS(1,ISING,XMESH(NMESH),UR,VR,ELAM)
CALL F06QFF('G',2,2,UR,2,URSTOR(1,1,NMESH),2)
CALL F06QFF('G',2,2,VR,2,VRSTOR(1,1,NMESH),2)
CALL F06QHF('G',2,2,0.D0,1.D0,RRSTOR(1,1,NMESH),2)
NRSTOR(NMESH) = 0.D0
ARGR = 0.D0
DO 60 I = NMESH,IMATCH + 1,-1
XO = XMESH(I)
XEND = XMESH(I-1)
CALL THTNT1(UR,VR,RINV,ARGR,NUX,PP(I),SP(I),QP(I),WP(I),ELAM,
& XO,XEND,PRN,IFAIL)
IF (IFAIL.NE.0) RETURN
DO 50 K = 1,2
DO 40 J = 1,2
URSTOR(J,K,I-1) = UR(J,K)
VRSTOR(J,K,I-1) = VR(J,K)
RRSTOR(J,K,I-1) = RINV(J,K)
40 CONTINUE
50 CONTINUE
NRSTOR(I-1) = NUX
60 CONTINUE
C Now match things up in the centre
CALL F06YAF('T','N',2,2,2,1.d0,VR,2,UL,2,0.D0,UC,2)
CALL F06YAF('T','N',2,2,2,-1.D0,UR,2,VL,2,1.D0,UC,2)
C Now determine the rank deficiency of UC and its left and
C right null-vectors.
MAXUC = 0.0D0
DO 80 I = 1,2
DO 70 J = 1,2
MAXUC = MAX(MAXUC,ABS(UC(I,J)))
70 CONTINUE
80 CONTINUE
IF (MAXUC.LT. (10.D0*X02AJF(0.D0))) THEN
MULTI = 2
DL1(1) = 1.D0
DL1(2) = 0.D0
DL2(1) = 0.D0
DL2(2) = 1.D0
DR1(1) = 1.D0
DR1(2) = 0.D0
DR2(1) = 0.D0
DR2(1) = 1.D0
ELSE IF ((UC(1,2)**2+UC(1,1)**2).GT. (UC(2,2)**2+UC(2,1)**2)) THEN
MULTI = 1
DL1(1) = -UC(1,2)
DL1(2) = UC(1,1)
ELSE
MULTI = 1
DL1(1) = -UC(2,2)
DL1(2) = UC(2,1)
END IF
IF (MAXUC.GT. (10.D0*X02AJF(0.D0))) THEN
IF ((UC(2,1)**2+UC(1,1)**2).GT. (UC(2,2)**2+UC(1,2)**2)) THEN
DR1(1) = -UC(2,1)
DR1(2) = UC(1,1)
ELSE
DR1(1) = -UC(2,2)
DR1(2) = UC(1,2)
END IF
END IF
C Now we have DL and DR we can calc Wi and Vi
LAST = .TRUE.
IF (MULTI.EQ.2) LAST = .FALSE.
90 IF (MULTI.EQ.2 .AND. LAST) THEN
C Set up conditions for second eigenfunction
WSTORE(1,IMATCH) = DL2(1)
WSTORE(2,IMATCH) = DL2(2)
ELSE
C Set up conditions for first eigenfunction
WSTORE(1,IMATCH) = DL1(1)
WSTORE(2,IMATCH) = DL1(2)
END IF
DO 100 I = IMATCH,1,-1
CALL F06PAF('N',2,2,1.D0,RLSTOR(1,1,I),2,WSTORE(1,I),1,0.D0,
& WSTORE(1,I-1),1)
C Scale to avoid possible accumulation to overflow:
DUMMY2 = SQRT(WSTORE(1,I-1)**2+WSTORE(2,I-1)**2)
IF (DUMMY2.GT.1.D0/SQRT(X02AJF(0.D0))) THEN
C DO 101 J = I-1,IMATCH
J = I - 1
WSTORE(1,J) = WSTORE(1,J)/DUMMY2
WSTORE(2,J) = WSTORE(2,J)/DUMMY2
NLSTOR(I) = NLSTOR(I) - LOG(DUMMY2)
C 101 CONTINUE
END IF
100 CONTINUE
IF (MULTI.EQ.2 .AND. LAST) THEN
C Set up conditions for second eigenfunction
VSTORE(1,IMATCH) = DR2(1)
VSTORE(2,IMATCH) = DR2(2)
ELSE
C Set up conditions for first eigenfunction
VSTORE(1,IMATCH) = DR1(1)
VSTORE(2,IMATCH) = DR1(2)
END IF
DO 110 I = IMATCH,NMESH - 1
CALL F06PAF('N',2,2,1.D0,RRSTOR(1,1,I),2,VSTORE(1,I),1,0.D0,
& VSTORE(1,I+1),1)
C Scale to avoid possible accumulation to overflow:
DUMMY2 = SQRT(VSTORE(1,I+1)**2+VSTORE(2,I+1)**2)
IF (DUMMY2.GT.1.D0/SQRT(X02AJF(0.D0))) THEN
C DO 111 J = IMATCH,I+1
J = I + 1
VSTORE(1,J) = VSTORE(1,J)/DUMMY2
VSTORE(2,J) = VSTORE(2,J)/DUMMY2
NRSTOR(J) = NRSTOR(J) - LOG(DUMMY2)
C 111 CONTINUE
END IF
110 CONTINUE
C Calculate Zl and Zr
CALL F06PAF('N',2,2,1.D0,ULSTOR(1,1,IMATCH),2,WSTORE(1,IMATCH),1,
& 0.D0,ZLSTOR(1,IMATCH),1)
CALL F06PAF('N',2,2,1.D0,VLSTOR(1,1,IMATCH),2,WSTORE(1,IMATCH),1,
& 0.D0,ZVLSTR(1,IMATCH),1)
EXCONT = 0.D0
DO 120 I = IMATCH - 1,0,-1
EXCONT = EXCONT + NLSTOR(I+1)
CALL F06PAF('N',2,2,EXP(-EXCONT),ULSTOR(1,1,I),2,WSTORE(1,I),
& 1,0.D0,ZLSTOR(1,I),1)
CALL F06PAF('N',2,2,EXP(-EXCONT),VLSTOR(1,1,I),2,WSTORE(1,I),
& 1,0.D0,ZVLSTR(1,I),1)
120 CONTINUE
CALL F06PAF('N',2,2,1.D0,URSTOR(1,1,IMATCH),2,VSTORE(1,IMATCH),1,
& 0.D0,ZRSTOR(1,IMATCH),1)
CALL F06PAF('N',2,2,1.D0,VRSTOR(1,1,IMATCH),2,VSTORE(1,IMATCH),1,
& 0.D0,ZVRSTR(1,IMATCH),1)
EXCONT = 0.D0
DO 130 I = IMATCH + 1,NMESH
EXCONT = EXCONT + NRSTOR(I-1)
CALL F06PAF('N',2,2,EXP(-EXCONT),URSTOR(1,1,I),2,VSTORE(1,I),
& 1,0.D0,ZRSTOR(1,I),1)
CALL F06PAF('N',2,2,EXP(-EXCONT),VRSTOR(1,1,I),2,VSTORE(1,I),
& 1,0.D0,ZVRSTR(1,I),1)
C WRITE (14,FMT=*)
C WRITE (14,FMT=*) 'URSTOR:'
C WRITE (14,FMT=*) URSTOR(1,1,I),URSTOR(1,2,I)
C WRITE (14,FMT=*) URSTOR(2,1,I),URSTOR(2,2,I)
C WRITE (14,FMT=*) 'VSTORE:'
C WRITE (14,FMT=*) VSTORE(1,I),VSTORE(2,I)
C WRITE (14,FMT=*)
C WRITE (14,FMT=*) 'I,EXCONT,NRSTOR(I-1),ZRSTOR(1,I)',I,EXCONT,
C & NRSTOR(I-1),ZRSTOR(1,I)
130 CONTINUE
C Now compute BL,BR
BETAL = 0.D0
DO 140 I = 1,IMATCH
BETAL = BETAL + WP(I)* (((ZLSTOR(1,I)+ZLSTOR(1,I-1))*0.5D0)**
& 2)* (XMESH(I)-XMESH(I-1))
140 CONTINUE
BETAR = 0.D0
DO 150 I = IMATCH + 1,NMESH
BETAR = BETAR + WP(I)* (((ZRSTOR(1,I)+ZRSTOR(1,I-1))*0.5D0)**
& 2)* (XMESH(I)-XMESH(I-1))
C WRITE (14,FMT=*) 'I,BETAR-so-far:',I,BETAR
C WRITE (14,FMT=*) 'WP(I),ZRSTOR(1,I),ZRSTOR(1,I-1),H:',
C & WP(I),ZRSTOR(1,I),ZRSTOR(1,I-1),XMESH(I)-XMESH(I-1)
150 CONTINUE
C Now calculat alfal and alfar
DUMMY(1) = ZLSTOR(1,IMATCH)**2 + ZRSTOR(1,IMATCH)**2
DUMMY(2) = ZLSTOR(2,IMATCH)**2 + ZRSTOR(2,IMATCH)**2
DUMMY(3) = ZVLSTR(1,IMATCH)**2 + ZVRSTR(1,IMATCH)**2
DUMMY(4) = ZVLSTR(2,IMATCH)**2 + ZVRSTR(2,IMATCH)**2
J = 1
DUMMY2 = DUMMY(J)
DO 160 JJ = 2,4
IF (DUMMY(JJ).GT.DUMMY2) THEN
J = JJ
DUMMY2 = DUMMY(JJ)
END IF
160 CONTINUE
IF (J.LE.2) THEN
DUMMY(1) = ZLSTOR(J,IMATCH)
DUMMY(2) = ZRSTOR(J,IMATCH)
ELSE
DUMMY(1) = ZVLSTR(J-2,IMATCH)
DUMMY(2) = ZVRSTR(J-2,IMATCH)
END IF
IF (DUMMY(1).GE.DUMMY(2)) THEN
RATIO = DUMMY(2)/DUMMY(1)
ALFAR = 1.D0/SQRT(BETAR+BETAL*RATIO**2)
ALFAL = ALFAR*RATIO
ELSE
RATIO = DUMMY(1)/DUMMY(2)
ALFAL = 1.D0/SQRT(BETAL+BETAR*RATIO**2)
ALFAR = ALFAL*RATIO
END IF
C WRITE (14,FMT=*) 'ALFAL,ALFAR:',ALFAL,ALFAR
C WRITE (14,FMT=*) 'DUMMY(1),DUMMY(2),DUMMY(3),DUMMY(4):',DUMMY(1),
C & DUMMY(2),DUMMY(3),DUMMY(4)
C WRITE (14,FMT=*) 'J:',J
C IF (J.LE.2) THEN
C WRITE (14,FMT=*) 'ZLSTOR(J,IMATCH),ZRSTOR(J,IMATCH):',
C & ZLSTOR(J,IMATCH),ZRSTOR(J,IMATCH)
C ELSE
C WRITE (14,FMT=*) 'ZVLSTR(J-2,IMATCH),ZVRSTR(J-2,IMATCH):',
C & ZVLSTR(J-2,IMATCH),ZVRSTR(J-2,IMATCH)
C END IF
C WRITE (14,FMT=*) 'BETAL,BETAR:',BETAL,BETAR
C Now compute yl(i) and yr(i)
IF (MULTI.EQ.2 .AND. LAST) THEN
DO 170 I = 0,IMATCH - 1
YSTOR2(1,I) = ALFAL*ZLSTOR(1,I)
YSTOR2(2,I) = ALFAL*ZLSTOR(2,I)
YSTOR2(3,I) = ALFAL*ZVLSTR(1,I)
YSTOR2(4,I) = ALFAL*ZVLSTR(2,I)
170 CONTINUE
DO 180 I = IMATCH,NMESH
YSTOR2(1,I) = ALFAR*ZRSTOR(1,I)
YSTOR2(2,I) = ALFAR*ZRSTOR(2,I)
YSTOR2(3,I) = ALFAR*ZVRSTR(1,I)
YSTOR2(4,I) = ALFAR*ZVRSTR(2,I)
180 CONTINUE
ELSE
DO 190 I = 0,IMATCH - 1
YSTOR1(1,I) = ALFAL*ZLSTOR(1,I)
YSTOR1(2,I) = ALFAL*ZLSTOR(2,I)
YSTOR1(3,I) = ALFAL*ZVLSTR(1,I)
YSTOR1(4,I) = ALFAL*ZVLSTR(2,I)
C WRITE (14,FMT=*) 'YSTOR1,ZLSTOR:',YSTOR1(1,I),ZLSTOR(1,I)
190 CONTINUE
DO 200 I = IMATCH,NMESH
YSTOR1(1,I) = ALFAR*ZRSTOR(1,I)
YSTOR1(2,I) = ALFAR*ZRSTOR(2,I)
YSTOR1(3,I) = ALFAR*ZVRSTR(1,I)
YSTOR1(4,I) = ALFAR*ZVRSTR(2,I)
C WRITE (14,FMT=*) 'YSTOR1,ZRSTOR:',YSTOR1(1,I),ZRSTOR(1,I)
200 CONTINUE
IF (MULTI.EQ.2) THEN
C LAST must be .FALSE., i.e. we have done only one
C eigenfunction.
LAST = .TRUE.
GO TO 90
END IF
END IF
RETURN
END
C ----------------------------------------------------------------------
SUBROUTINE THTNT1(U,V,RINV,ARGDET,NUX,P,S,Q,W,ELAM,XO,XEND,PRN,
& IFAIL)
C This subroutine integrates the theta matrix of a 4th order
C Sturm-Liouville problem across an interval on which the
C coefficients are constant, updating argdet(theta) in the
C process. The theta matrix itself is represented by two
C 2x2 matrices u and v such that theta = (v+iu)(v-iu)^{-1}.
C First we update the theta matrix (i.e. the matrices u and
C v) so that we will have all the necessary information for
C updating argdet(theta).
C The formulae which we use for doing the update are as follows.
C If xo < xend, then
C
C Argdet(Theta(xend)) = Argdet(Theta(xo))
C + (Sum of Phase Angles of Theta(xend) in [0,2pi))
C - (Sum of Phase Angles of Theta_{0}^{*}Theta(xo) in [0,2pi))
C - (Sum of Phase Angles of Theta_{0}(xo) in [0,2pi))
C + 2(n+2)pi,
C
C where n is the number of zeros of det(Uo) in [xo,xend) (remember
C that Theta_{0}(xend) = I).
C If xo > xend, then
C
C Argdet(Theta(xend)) = Argdet(Theta(xo))
C - (Sum of Phase Angles of Theta^{*}(xend) in [0,2pi))
C + (Sum of Phase Angles of Theta^{*}Theta_{0}(xo) in [0,2pi))
C - (Sum of Phase Angles of Theta_{0}(xo) in (0,2pi])
C - 2.n.pi
C
C .. Scalar Arguments ..
DOUBLE PRECISION ARGDET,ELAM,NUX,P,Q,S,W,XEND,XO
INTEGER IFAIL
LOGICAL PRN
C ..
C .. Array Arguments ..
DOUBLE PRECISION RINV(2,2),U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,ARG0,ARGEND,ARGTHO,B,C,D,DISC1,DISC2,DUMMY,EC,
& M,OMEGA1,OMEGA2,SCAL,TWOPI
INTEGER I,J,N
LOGICAL LINTSP,LNEG,LPOS,SPCIAL
C ..
C .. Local Arrays ..
DOUBLE PRECISION CU(2,2),CV(2,2),RNORM(2,2),SU(2,2),SV(2,2),
& ULOC(2,2),UO(2,2),UTEMP(2,2),VLOC(2,2),VO(2,2),
& VTEMP(2,2),Y1(4),Y2(4),Y3(4),Y4(4)
C ..
C .. External Functions ..
DOUBLE PRECISION X01AAF
EXTERNAL X01AAF
C ..
C .. External Subroutines ..
EXTERNAL F06YAF,FUNSOL,ITNLM1,LRGLM1,NEGLM1,SL4CNT,SPTH4
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,DBLE,LOG,MAX,SIN,SQRT
C ..
TWOPI = 2.D0*X01AAF(0.D0)
C Copy the matrices into temporary storage for later
C use.
DO 20 J = 1,2
DO 10 I = 1,2
UTEMP(I,J) = U(I,J)
VTEMP(I,J) = V(I,J)
10 CONTINUE
20 CONTINUE
C
DISC1 = 4.D0*P* (ELAM*W-Q)
DISC2 = S*S + DISC1
LPOS = DISC1 .GE. 0.D0
LINTSP = DISC2 .GE. 0.D0 .AND. S - SQRT(ABS(DISC2)) .GE. 0.D0
LNEG = DISC2 .LT. 0.D0
IF (LPOS .OR. LINTSP) THEN
SCAL = ABS(XEND-XO)*SQRT(0.5D0* (S+SQRT(ABS(DISC2)))/P)
ELSE IF (LNEG) THEN
SCAL = ABS(XEND-XO)*0.5D0*SQRT((S+SQRT(ABS(DISC1)))/P)
ELSE
SCAL = 0.D0
END IF
SPCIAL = SCAL .GT. 1.D-1
IF (.NOT.SPCIAL) THEN
C In this case we must accept the output from FUNSOL as correct
C WRITE (8,FMT=*) 'Doing simple case'
CALL FUNSOL(XO,XEND,P,S,Q,W,ELAM,Y1,Y2,Y3,Y4,SCAL)
CU(1,1) = Y1(1)
CU(1,2) = Y2(1)
CU(2,1) = Y1(2)
CU(2,2) = Y2(2)
SU(1,1) = Y3(1)
SU(1,2) = Y4(1)
SU(2,1) = Y3(2)
SU(2,2) = Y4(2)
SV(1,1) = Y1(3)
SV(1,2) = Y2(3)
SV(2,1) = Y1(4)
SV(2,2) = Y2(4)
CV(1,1) = Y3(3)
CV(1,2) = Y4(3)
CV(2,1) = Y3(4)
CV(2,2) = Y4(4)
C
C Now do the update of Theta from xo to xend, using the formulae
C u(xend) = cu.u(xo)+su.v(xo),
C v(xend) = sv.u(xo)+cv.v(xo).
C
CALL F06YAF('N','N',2,2,2,1.D0,CU,2,UTEMP,2,0.D0,UO,2)
CALL F06YAF('N','N',2,2,2,1.D0,SU,2,VTEMP,2,1.D0,UO,2)
CALL F06YAF('N','N',2,2,2,1.D0,CV,2,VTEMP,2,0.D0,VO,2)
CALL F06YAF('N','N',2,2,2,1.D0,SV,2,UTEMP,2,1.D0,VO,2)
C
C Now do some postmultiplication to avoid ill-conditioning
C of these matrices.
C
B = -UO(1,1)*VO(2,2) + UO(1,2)*VO(2,1) - UO(2,2)*VO(1,1) +
& UO(2,1)*VO(1,2)
AC = UO(1,1)*UO(2,2) - UO(2,1)*UO(1,2) - VO(1,1)*VO(2,2) +
& VO(2,1)*VO(1,2)
C = UO(1,1)*UO(2,2) - UO(2,1)*UO(1,2) + VO(1,1)*VO(2,2) -
& VO(2,1)*VO(1,2)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C = COS(M/2.D0)
AC = SIN(M/2.D0)
C Form the normalising matrix which is equal to the return velaue of R:
DO 40 J = 1,2
DO 30 I = 1,2
RNORM(I,J) = UO(I,J)*C - VO(I,J)*AC
30 CONTINUE
40 CONTINUE
C In this case a uniform scaling is used by the routine FUNSOL, so
NUX = SCAL
C Form the inverse of the normalising matrix:
D = RNORM(1,1)*RNORM(2,2) - RNORM(1,2)*RNORM(2,1)
EC = RNORM(1,1)/D
RINV(1,1) = RNORM(2,2)/D
RINV(2,2) = EC
RINV(1,2) = -RNORM(1,2)/D
RINV(2,1) = -RNORM(2,1)/D
C
C Peform the postmultiplication:
CALL F06YAF('N','N',2,2,2,1.D0,UO,2,RINV,2,0.D0,U,2)
CALL F06YAF('N','N',2,2,2,1.D0,VO,2,RINV,2,0.D0,V,2)
C U and V have now been postmultiplied and V+iU should be
C reasonably well-conditioned.
C Now set the U and V for Theta_{0}
UO(1,1) = -SU(1,1)
UO(1,2) = SU(1,2)
UO(2,1) = SU(2,1)
UO(2,2) = -SU(2,2)
VO(1,1) = CV(1,1)
VO(1,2) = -CV(1,2)
VO(2,1) = -CV(2,1)
VO(2,2) = CV(2,2)
C WRITE (8,FMT=*) XO,XEND,'Used FUNSOL'
ELSE
IF (LPOS) THEN
CALL LRGLM1(ELAM,NUX,RINV,P,S,Q,W,U,V,UO,VO,XO,XEND)
C WRITE (8,FMT=*) XO,XEND,'Used LRGLM1'
ELSE IF (LINTSP) THEN
CALL ITNLM1(ELAM,NUX,RINV,P,S,Q,W,U,V,UO,VO,XO,XEND)
C WRITE (8,FMT=*) XO,XEND,'Used ITNLM1'
ELSE IF (LNEG) THEN
CALL NEGLM1(ELAM,NUX,RINV,P,S,Q,W,U,V,UO,VO,XO,XEND)
C WRITE (8,FMT=*) XO,XEND,'Used NEGLM1'
END IF
END IF
C Normalise U and V so that the greatest element of the two is
C equal to 1.
DUMMY = 0.D0
DO 60 J = 1,2
DO 50 I = 1,2
DUMMY = MAX(DUMMY,ABS(U(I,J)),ABS(V(I,J)))
50 CONTINUE
60 CONTINUE
DO 80 J = 1,2
DO 70 I = 1,2
U(I,J) = U(I,J)/DUMMY
V(I,J) = V(I,J)/DUMMY
70 CONTINUE
80 CONTINUE
C Include this into the value of NUX
NUX = NUX + LOG(DUMMY)
C Now all we have to do is get the eigenvalues of lots of
C different matrices.
C Get the eigenvalues of Theta_{0}^{*}Theta(xo).
CALL F06YAF('T','N',2,2,2,1.D0,VO,2,UTEMP,2,0.D0,ULOC,2)
CALL F06YAF('T','N',2,2,2,-1.D0,UO,2,VTEMP,2,1.D0,ULOC,2)
CALL F06YAF('T','N',2,2,2,1.D0,VO,2,VTEMP,2,0.D0,VLOC,2)
CALL F06YAF('T','N',2,2,2,1.D0,UO,2,UTEMP,2,1.D0,VLOC,2)
CALL SPTH4(ULOC,VLOC,OMEGA1,OMEGA2,ARG0,'L',IFAIL)
IF (IFAIL.NE.0) THEN
IFAIL = 5
RETURN
END IF
IF (XO.GT.XEND) THEN
C We actually want the phase-angles of the hermitian conjugate
C of Theta_{0}^{*}Theta(xo)
IF (OMEGA1.NE.0.D0) OMEGA1 = TWOPI - OMEGA1
IF (OMEGA2.NE.0.D0) OMEGA2 = TWOPI - OMEGA2
ARG0 = -OMEGA1 - OMEGA2
END IF
C Get the phase angles of Theta(xend)
CALL SPTH4(U,V,OMEGA1,OMEGA2,ARGEND,'L',IFAIL)
IF (IFAIL.NE.0) THEN
IFAIL = 6
RETURN
END IF
IF (XO.GT.XEND) THEN
C We actually want the phase-angles of the hermitian conjugate
C of Theta(xend)
IF (OMEGA1.NE.0.D0) OMEGA1 = TWOPI - OMEGA1
IF (OMEGA2.NE.0.D0) OMEGA2 = TWOPI - OMEGA2
ARGEND = -OMEGA1 - OMEGA2
END IF
C Get phase-angles of Theta_{0}(xo), normalised to lie in (0,2pi]
CALL SPTH4(UO,VO,OMEGA1,OMEGA2,ARGTHO,'R',IFAIL)
IF (IFAIL.NE.0) THEN
IFAIL = 7
RETURN
END IF
C Do zero count:
CALL SL4CNT(XEND-XO,P,S,Q,W,ELAM,N)
C Do the update of ARGDET:
IF (XO.LE.XEND) THEN
ARGDET = ARGDET + ARGEND - ARG0 + DBLE(N+2)*TWOPI - ARGTHO
ELSE
ARGDET = ARGDET + ARGEND - ARG0 - DBLE(N)*TWOPI - ARGTHO
END IF
END
C --------------------------------------------------------------------
SUBROUTINE LRGLM1(ELAM,NUX,RINV,P,S,Q,W,U,V,SU,CV,XO,XEND)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,NUX,P,Q,S,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION CV(2,2),RINV(2,2),SU(2,2),U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,AOMX,B,BETA,C,COSHCS,COSHNX,COSHSN,COSOMX,
& DETCU,DETCV,DETR,DETSU,DETSV,DETU,DETV,DISC,FAC,
& H,H2,H2SNSN,M,NO,NO2,NU,NU2,NU3,NU4,NU6,OMEGA,
& OMEGA2,OMEGA3,OMEGA4,OMEGA6,OMX,PNO,PNO2,PNORAD,
& PRAD,PRAD2,RAD,RAD2,SINHCS,SINHNX,SINHSN,SINOMX,
& TRUVA,ZCOSMT,ZSINMT
INTEGER I,J
C ..
C .. Local Arrays ..
DOUBLE PRECISION A(2,2),CU(2,2),CUSVA(2,2),R(2,2),SUCVA(2,2),
& SV(2,2),T(2,2),TEMP(2,2),UT(2,2),UUA(2,2),
& UVA(2,2),VT(2,2),VUA(2,2),VVA(2,2)
C ..
C .. External Functions ..
DOUBLE PRECISION PHI
EXTERNAL PHI
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,EXP,SIN,SQRT
C ..
C .. External Subroutines ..
EXTERNAL F06YAF
C ..
DISC = SQRT(S*S+4.D0*P* (ELAM*W-Q))
NU = SQRT(0.5d0* (S+DISC)/P)
OMEGA = SQRT(0.5d0* (-S+DISC)/P)
DETU = U(1,1)*U(2,2) - U(2,1)*U(1,2)
DETV = V(1,1)*V(2,2) - V(2,1)*V(1,2)
H = XEND - XO
H2 = H*H
NUX = ABS(NU*H)
FAC = EXP(-NUX)
IF (NUX.GE.0.1D0) THEN
SINHNX = (1.D0-FAC*FAC)/ (2.D0*NUX)
ELSE
SINHNX = PHI(NUX*NUX)*FAC
END IF
COSHNX = (1.D0+FAC*FAC)/2.D0
OMX = OMEGA*H
AOMX = ABS(OMX)
IF (AOMX.GE.0.1D0) THEN
SINOMX = SIN(OMX)/OMX
ELSE
SINOMX = PHI(-OMX*OMX)
END IF
COSOMX = COS(OMX)
COSHSN = H*COSHNX*SINOMX
SINHCS = H*SINHNX*COSOMX
COSHCS = COSHNX*COSOMX
SINHSN = SINHNX*SINOMX
H2SNSN = H2*SINHSN
OMEGA2 = OMEGA*OMEGA
NU2 = NU*NU
OMEGA3 = OMEGA*OMEGA2
NU3 = NU*NU2
OMEGA4 = OMEGA2*OMEGA2
NU4 = NU2*NU2
OMEGA6 = OMEGA2*OMEGA4
NU6 = NU2*NU4
RAD = OMEGA2 + NU2
RAD2 = RAD*RAD
PRAD = P*RAD
PRAD2 = PRAD*PRAD
NO = NU*OMEGA
PNO = P*NO
PNORAD = PNO*RAD
NO2 = NO*NO
PNO2 = P*NO2
C Compute the CU,CV,SU,SV matrices here
CU(1,1) = (NU2*COSOMX*FAC+OMEGA2*COSHNX)/RAD
CU(1,2) = H* (NU2*SINHNX+OMEGA2*SINOMX*FAC)/RAD
CU(2,1) = H* (-NU2*OMEGA2*SINOMX*FAC+OMEGA2*NU2*SINHNX)/RAD
CU(2,2) = (NU2*COSHNX+OMEGA2*COSOMX*FAC)/RAD
SU(1,1) = H* (SINOMX*FAC-SINHNX)/PRAD
SU(1,2) = (COSHNX-COSOMX*FAC)/PRAD
SU(2,1) = -SU(1,2)
SU(2,2) = H* (NU2*SINHNX+OMEGA2*SINOMX*FAC)/PRAD
CV(1,1) = (NU2*COSOMX*FAC+OMEGA2*COSHNX)/RAD
CV(1,2) = H* (NO* (NO*SINOMX*FAC-NO*SINHNX))/RAD
CV(2,1) = H* (-OMEGA2*SINOMX*FAC-NU2*SINHNX)/RAD
CV(2,2) = (NU2*COSHNX+OMEGA2*COSOMX*FAC)/RAD
SV(1,1) = H* (-PNO* (NU2*NO*SINOMX*FAC+OMEGA2*NO*SINHNX))/RAD
SV(1,2) = (-P*NU2*OMEGA2* (COSHNX-COSOMX*FAC))/RAD
SV(2,1) = -SV(1,2)
SV(2,2) = H* (P* (NU4*SINHNX-OMEGA4*SINOMX*FAC))/RAD
C Now we can compute UT and VT
CALL F06YAF('N','N',2,2,2,1.D0,CU,2,U,2,0.D0,UT,2)
CALL F06YAF('N','N',2,2,2,1.D0,SU,2,V,2,1.D0,UT,2)
CALL F06YAF('N','N',2,2,2,1.D0,SV,2,U,2,0.D0,VT,2)
CALL F06YAF('N','N',2,2,2,1.D0,CV,2,V,2,1.D0,VT,2)
C Compute SU.Adjoint(CV).exp(-NUX):
SUCVA(1,1) = (COSHSN-SINHCS)/PRAD
SUCVA(1,2) = ((OMEGA2-NU2)* (FAC-COSHCS)+2.D0*NO2*H2SNSN)/
& (PRAD*RAD)
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = (NU2*SINHCS+OMEGA2*COSHSN)/PRAD
C Compute CU.Adjoint(SV).exp(-NUX):
CUSVA(1,1) = P* (NU4*SINHCS-OMEGA4*COSHSN)/RAD
CUSVA(1,2) = -PNO2* ((NU4+OMEGA4)*H2SNSN+
& (NU2-OMEGA2)* (FAC-COSHCS))/RAD2
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = -PNO2* (NU2*COSHSN+OMEGA2*SINHCS)/RAD
C Compute A = U.Adjoint(V):
A(1,1) = U(1,1)*V(2,2) - U(1,2)*V(2,1)
A(1,2) = U(1,2)*V(1,1) - U(1,1)*V(1,2)
A(2,1) = A(1,2)
A(2,2) = V(1,1)*U(2,2) - U(2,1)*V(1,2)
C Compute T = (CU.A.Adjoint(CV) + SU.A.Adjoint(SV))*exp(-NUX):
T(1,1) = A(1,1)*COSHCS + 2.D0*A(1,2)* (OMEGA2*COSHSN+NU2*SINHCS)/
& RAD + A(2,2)*H2SNSN
T(1,2) = (A(1,1)*NO2* (SINHCS-COSHSN)+
& A(1,2)* ((OMEGA2-NU2)* (FAC* (OMEGA2-NU2)+
& 2.D0*NO2*H2SNSN)+4.D0*NO2*COSHCS)/RAD+
& A(2,2)* (NU2*SINHCS+OMEGA2*COSHSN))/RAD
T(2,1) = T(1,2)
T(2,2) = -A(1,1)*NO2*H2SNSN + 2.D0*A(1,2)*NO2* (SINHCS-COSHSN)/
& RAD + A(2,2)*COSHCS
C Compute U.Adjoint(V) = (CU.Uo+SU.Vo).Adjoint(SV.Uo+CV.Vo):
UVA(1,1) = DETU*CUSVA(1,1) + DETV*SUCVA(1,1) + T(1,1)
UVA(1,2) = DETU*CUSVA(1,2) + DETV*SUCVA(1,2) + T(1,2)
UVA(2,1) = DETU*CUSVA(2,1) + DETV*SUCVA(2,1) + T(2,1)
UVA(2,2) = DETU*CUSVA(2,2) + DETV*SUCVA(2,2) + T(2,2)
VUA(1,1) = UVA(2,2)
VUA(2,2) = UVA(1,1)
VUA(1,2) = -UVA(1,2)
VUA(2,1) = -UVA(2,1)
DETCU = ((NU4+OMEGA4)*COSHCS+NO2* ((NU2-OMEGA2)*H2SNSN+2.D0*FAC))/
& RAD2
DETSU = ((NU2-OMEGA2)*H2SNSN+2.D0* (FAC-COSHCS))/PRAD2
DETCV = DETCU
DETSV = P*PNO* ((OMEGA6-NU6)*NO*H2SNSN+
& 2.D0*NU3*OMEGA3* (FAC-COSHCS))/RAD2
C Compute trace(CU.A.Adjoint(SU)).exp(-NUX):
BETA = A(1,1)* (OMEGA2*COSHSN+NU2*SINHCS)/PRAD -
& 2.D0*A(1,2)* ((NU2-OMEGA2)* (FAC-COSHCS)-2.D0*NO2*H2SNSN)/
& (P*RAD2) + A(2,2)* (COSHSN-SINHCS)/PRAD
UUA(1,1) = DETU*DETCU + DETV*DETSU + BETA
UUA(1,2) = 0.D0
UUA(2,1) = 0.D0
UUA(2,2) = UUA(1,1)
C Compute trace(CV.A.Adjoint(SV)).exp(-NUX):
BETA = (-A(1,1)*PNO2* (NU2*COSHSN+OMEGA2*SINHCS)-
& 2.D0*PNO2*A(1,2)* ((NU2-OMEGA2)* (FAC-COSHCS)+ (NU4+
& OMEGA4)*H2SNSN)/RAD-A(2,2)*P* (OMEGA4*COSHSN-NU4*SINHCS))/
& RAD
VVA(1,1) = DETU*DETSV + DETV*DETCV + BETA
VVA(1,2) = 0.D0
VVA(2,1) = 0.D0
VVA(2,2) = VVA(1,1)
TRUVA = UVA(1,1) + UVA(2,2)
AC = UUA(1,1) - VVA(1,1)
B = -TRUVA
C = UUA(1,1) + VVA(1,1)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C Now we have got the value of M we can compute R.
ZCOSMT = COS(M/2.D0)
ZSINMT = SIN(M/2.D0)
DO 20 J = 1,2
DO 10 I = 1,2
R(I,J) = UT(I,J)*ZCOSMT - VT(I,J)*ZSINMT
10 CONTINUE
20 CONTINUE
C Compute the DET R and RINV
DETR = ZCOSMT**2*UUA(1,1) + ZSINMT**2*VVA(1,1) -
& ZSINMT*ZCOSMT* (UVA(1,1)+UVA(2,2))
RINV(1,1) = R(2,2)/DETR
RINV(1,2) = -R(1,2)/DETR
RINV(2,1) = -R(2,1)/DETR
RINV(2,2) = R(1,1)/DETR
CALL F06YAF('N','N',2,2,2,1.D0,UT,2,RINV,2,0.D0,TEMP,2)
DO 40 J = 1,2
DO 30 I = 1,2
U(I,J) = (UUA(I,J)*ZCOSMT-UVA(I,J)*ZSINMT)/DETR
V(I,J) = (VUA(I,J)*ZCOSMT-VVA(I,J)*ZSINMT)/DETR
30 CONTINUE
40 CONTINUE
NUX = 0.D0
C That's done the very tricky renormalisation. Now do the SU and CV
C matrices which are the U and V matrices for Theta_{0}.
TRUVA = SUCVA(1,1) + SUCVA(2,2)
AC = (DETSU-DETCV)
B = -TRUVA
C = (DETSU+DETCV)
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C NOTE: in the following formulae we want CV and SU for BACKWARDS
C integration. We exploit the fact that the diagonal terms of
C SUCVA are odd functions of the direction of integration while
C the off-diagonal terms are even functions of the direction
C of integration, while the determinants which appear in the
C formulae are even functions of the direction of integration.
SU(1,1) = DETSU*COS(M/2.D0) + SUCVA(1,1)*SIN(M/2.D0)
SU(1,2) = -SUCVA(1,2)*SIN(M/2.D0)
SU(2,1) = -SUCVA(2,1)*SIN(M/2.D0)
SU(2,2) = DETSU*COS(M/2.D0) + SUCVA(2,2)*SIN(M/2.D0)
CV(1,1) = -SUCVA(2,2)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
CV(1,2) = -SUCVA(1,2)*COS(M/2.D0)
CV(2,1) = -SUCVA(2,1)*COS(M/2.D0)
CV(2,2) = -SUCVA(1,1)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE ITNLM1(ELAM,NUX,RINV,P,S,Q,W,U,V,SU,CV,XO,XEND)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,NUX,P,Q,S,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION CV(2,2),RINV(2,2),SU(2,2),U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,ADJUST,AEHH,AETAH,ARHH,B,BETA,C,CEHCHH,CEHSHH,
& CMSN,CNCM,CNSM,COSHEH,COSHHH,COSHMX,COSHNX,DETCU,
& DETCV,DETR,DETSU,DETSV,DETU,DETV,DISC,ETA,ETA2,
& ETAH,ETASQ,ETSMRS,ETSRHS,FAC,FACL,H,H2,M,MU,MU2,
& MU3,MU4,MU6,MUX,NM,NU,NU2,NU3,NU4,NU6,PNM,PNMR,
& PRAD,PRAD2,RAD,RAD2,RADH2,RH,RH2,RHH,RHSQ,SEHCHH,
& SEHSHH,SINHMX,SINHNX,SN2ERT,SN2HRT,SNERAT,SNERNS,
& SNHR1,SNHR2,SNHRAT,SNHRNS,SNSM,TRUVA,ZCOSMT,
& ZSINMT
INTEGER I,J
LOGICAL SMALL
C ..
C .. Local Arrays ..
C DOUBLE PRECISION CUT(2,2),CVT(2,2),SUT(2,2),SVT(2,2)
DOUBLE PRECISION A(2,2),CU(2,2),CUSVA(2,2),R(2,2),SUCVA(2,2),
& SV(2,2),T(2,2),UT(2,2),UUA(2,2),UVA(2,2),VT(2,2),
& VUA(2,2),VVA(2,2)
C ..
C .. External Functions ..
DOUBLE PRECISION PHI,SNHDIF
EXTERNAL PHI,SNHDIF
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,EXP,MIN,SIGN,SIN,SQRT
C ..
C .. External Subroutines ..
EXTERNAL F06YAF
C ..
C WRITE (8,FMT=*) 'Debugging output from ITNLM1'
DISC = SQRT(S*S+4.D0*P* (ELAM*W-Q))
NU2 = 0.5D0* (S+DISC)/P
MU2 = NU2 - DISC/P
NU = SQRT(NU2)
MU = SQRT(MU2)
NU3 = NU*NU2
MU3 = MU*MU2
NU4 = NU2*NU2
MU4 = MU2*MU2
NU6 = NU3*NU3
MU6 = MU3*MU3
DETU = U(1,1)*U(2,2) - U(2,1)*U(1,2)
DETV = V(1,1)*V(2,2) - V(2,1)*V(1,2)
H = XEND - XO
H2 = H*H
NUX = ABS(NU*H)
SINHNX = (1.D0-EXP(-2.D0*NUX))/2.D0
COSHNX = 1.D0 - SINHNX
SINHNX = SINHNX*SIGN(1.D0,H)
MUX = ABS(MU*H)
SINHMX = (1.D0-EXP(-2.D0*MUX))/2.D0
COSHMX = 1.D0 - SINHMX
SINHMX = SINHMX*SIGN(1.D0,H)
FAC = EXP(-NUX-MUX)
CMSN = COSHMX*SINHNX
CNSM = COSHNX*SINHMX
CNCM = COSHNX*COSHMX
SNSM = SINHNX*SINHMX
C If any of the following quantities is small then we are
C in trouble; Taylor series expansions have to be used.
C RAD = NU**2-MU**2
RAD = DISC/P
RAD2 = RAD*RAD
PRAD = DISC
PRAD2 = DISC*DISC
NM = NU*MU
PNM = P*NM
PNMR = PNM*RAD
SMALL = H2*MIN(ABS(RAD),ABS(PRAD),ABS(PNMR),PNM) .LT. 1.D-6
C WRITE (8,FMT=*) 'SMALL is ',SMALL
C Compute SU.Adjoint(CV).exp(-(NUX+MUX)):
IF (SMALL) THEN
C Compute certain quantities which will be jolly useful
C when calculating the matrix elements
ETA2 = NU + MU
ETA = 0.5D0*ETA2
ETASQ = ETA*ETA
ETAH = ETA*H
AETAH = ABS(ETAH)
RH2 = NU - MU
RH = 0.5D0*RH2
RHSQ = RH*RH
RHH = RH*H
ARHH = ABS(RHH)
COSHEH = 0.5D0* (1.D0+EXP(-2.D0*AETAH))
COSHHH = 0.5D0* (1.D0+EXP(-2.D0*ARHH))
CEHCHH = COSHEH*COSHHH
IF (ARHH.LT.0.1D0) THEN
SNHRAT = PHI(RHH*RHH)
SNHRNS = SNHRAT*EXP(-ARHH)
SNHRAT = SNHRAT*SNHRAT*FAC
SN2HRT = PHI(4.D0*RHH*RHH)*FAC
ELSE
FACL = EXP(-MUX)
SNHRAT = 0.5D0* (1.D0-EXP(-2.D0*ARHH))/ARHH
SNHRNS = SNHRAT
SNHRAT = SNHRAT*FACL
SNHRAT = SNHRAT*SNHRAT
SN2HRT = (1.D0-EXP(-4.D0*ARHH))/ (4.D0*ARHH)
SN2HRT = SN2HRT*FACL*FACL
END IF
IF (AETAH.LT.0.1D0) THEN
SNERAT = PHI(ETAH*ETAH)
SNERNS = SNERAT*EXP(-AETAH)
SNERAT = SNERAT*SNERAT*FAC
SN2ERT = PHI(4.D0*ETAH*ETAH)*FAC
ELSE
SNERAT = 0.5D0* (1.D0-EXP(-2.D0*AETAH))/AETAH
SNERNS = SNERAT
SNERAT = SNERAT*SNERAT
SN2ERT = (1.D0-EXP(-4.D0*AETAH))/ (4.D0*AETAH)
END IF
C WRITE (8,FMT=*) 'Check COSHEH:',COSH(ETAH)*EXP(-AETAH)-COSHEH
C WRITE (8,FMT=*) 'Check COSHRH:',COSH(RHH)*EXP(-ARHH)-COSHHH
C WRITE (8,FMT=*) 'Check SNHRNS:',SINH(RHH)*EXP(-ARHH)-
C & SNHRNS*RHH
C WRITE (8,FMT=*) 'Check SNERNS:',SINH(ETAH)*EXP(-ARHH)-
C & SNERNS*ETAH
CEHSHH = COSHEH*SNHRNS*H
SEHSHH = SNERNS*SNHRNS*H2
SEHCHH = COSHHH*SNERNS*H
IF (MUX.LT.0.1D0) THEN
SNHR1 = PHI(MUX*MUX)*EXP(-MUX)
ELSE
SNHR1 = SIGN(1.D0,H)*SINHMX/MUX
END IF
IF (NUX.LT.0.1D0) THEN
SNHR2 = PHI(NUX*NUX)*EXP(-NUX)
ELSE
SNHR2 = SIGN(1.D0,H)*SINHNX/NUX
END IF
END IF
C Here we can compute the SU,SV,CU,CV matrices
ADJUST = EXP(-NUX+MUX)
IF (.NOT.SMALL) THEN
CU(1,1) = (NU2*COSHMX*ADJUST-MU2*COSHNX)/RAD
CU(1,2) = (NU*SINHNX-MU*SINHMX*ADJUST)/RAD
CU(2,1) = NM* (NU*SINHMX*ADJUST-MU*SINHNX)/RAD
CU(2,2) = (NU2*COSHNX-MU2*COSHMX*ADJUST)/RAD
SU(1,1) = (((SINHMX)/MU)*ADJUST- ((SINHNX)/NU))/PRAD
SU(1,2) = (COSHNX-COSHMX*ADJUST)/PRAD
SU(2,1) = -SU(1,2)
SU(2,2) = (NU*SINHNX-MU*SINHMX*ADJUST)/PRAD
SV(1,1) = PNM* (NU3*SINHMX*ADJUST-MU3*SINHNX)/RAD
SV(1,2) = PNM*NM* (COSHNX-COSHMX*ADJUST)/RAD
SV(2,1) = -SV(1,2)
SV(2,2) = (P* (NU3*SINHNX-MU3*SINHMX*ADJUST))/RAD
CV(1,1) = (NU2*COSHMX*ADJUST-MU2*COSHNX)/RAD
CV(1,2) = -NM* (NU*SINHMX*ADJUST-MU*SINHNX)/RAD
CV(2,1) = - (NU*SINHNX-MU*SINHMX*ADJUST)/RAD
CV(2,2) = (NU2*COSHNX-MU2*COSHMX*ADJUST)/RAD
ELSE
ETSRHS = ETASQ + RHSQ
ETSMRS = ETASQ - RHSQ
CU(1,1) = CEHCHH - 0.5D0*ETSRHS*SEHSHH
CU(1,2) = 0.5D0* (CEHSHH+SEHCHH)
CU(2,1) = 0.5D0*ETSMRS* (SEHCHH-CEHSHH)
CU(2,2) = CEHCHH + 0.5D0*ETSRHS*SEHSHH
SU(1,2) = (0.5D0*SEHSHH)/P
SU(2,1) = -SU(1,2)
SU(2,2) = (0.5D0* (CEHSHH+SEHCHH))/P
CV(1,1) = CU(1,1)
CV(1,2) = -CU(2,1)
CV(2,1) = -CU(1,2)
CV(2,2) = CU(2,2)
SV(1,1) = 0.5D0*PNM* ((3.D0*ETASQ+RHSQ)*SEHCHH-
& (ETASQ+3.D0*RHSQ)*CEHSHH)
SV(1,2) = 0.5D0*PNM*NM*SEHSHH
SV(2,1) = -SV(1,2)
SV(2,2) = 0.5D0*P* ((ETASQ+3.0D0*RHSQ)*CEHSHH+
& (RHSQ+3.0D0*ETASQ)*SEHCHH)
RADH2 = ABS(RAD)*H2
AEHH = ABS(ETASQ-RHSQ)*H2
IF (RADH2.GE.0.01D0 .AND. AEHH.LT.0.01D0) THEN
C WRITE (8,FMT=*) 'Case 1 computation of SU(1,1)'
SU(1,1) = (1/PRAD)* (SNHR1*ADJUST-SNHR2)*H
ELSE IF (RADH2.LT.0.01D0 .AND. AEHH.GE.0.01D0) THEN
C WRITE (8,FMT=*) 'Case 2 computation of SU(1,1)'
SU(1,1) = H* (SEHCHH-CEHSHH)/ (2.D0*PNM)
ELSE
C WRITE (8,FMT=*) 'Case 3 computation of SU(1,1)'
SU(1,1) = ((H**3)/P)*SNHDIF(NUX,MUX)*EXP(-NUX)
END IF
END IF
C Now we can compute UT and VT
CALL F06YAF('N','N',2,2,2,1.D0,CU,2,U,2,0.D0,UT,2)
CALL F06YAF('N','N',2,2,2,1.D0,SU,2,V,2,1.D0,UT,2)
CALL F06YAF('N','N',2,2,2,1.D0,SV,2,U,2,0.D0,VT,2)
CALL F06YAF('N','N',2,2,2,1.D0,CV,2,V,2,1.D0,VT,2)
IF (.NOT.SMALL) THEN
SUCVA(1,1) = (NU*CNSM-MU*CMSN)/PNMR
SUCVA(1,2) = - (S* (FAC-CNCM)+2.D0*PNM*SNSM)/PRAD2
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = (NU*CMSN-MU*CNSM)/PRAD
ELSE
RADH2 = ABS(RAD)*H2
AEHH = ABS(ETASQ-RHSQ)*H2
IF (RADH2.GE.0.01D0 .AND. AEHH.LT.0.01D0) THEN
SUCVA(1,1) = H* (COSHNX*SNHR1-COSHMX*SNHR2)/PRAD
ELSE IF (RADH2.LT.0.01D0 .AND. AEHH.GE.0.01D0) THEN
SUCVA(1,1) = H* (SN2ERT-SN2HRT)/ (2.D0*PNM)
ELSE
SUCVA(1,1) = 2.D0* (H**3)*FAC*
& SNHDIF(2.D0*AETAH,2.D0*ARHH)/P
END IF
SUCVA(1,2) = 0.25D0*H*H* (SNERAT+SNHRAT)/P
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = 0.5D0*H* (SN2ERT+SN2HRT)/P
END IF
C Compute CU.Adjoint(SV).exp(-NUX-MUX):
IF (.NOT.SMALL) THEN
CUSVA(1,1) = P* (NU3*CMSN-MU3*CNSM)/RAD
CUSVA(1,2) = PNM* (NM* (NU2+MU2)* (FAC-CNCM)+ (NU4+MU4)*SNSM)/
& RAD2
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = PNM* (NU3*CNSM-MU3*CMSN)/RAD
ELSE
FACL = 0.25D0*PNM*H
CUSVA(1,1) = P*H*0.5D0* ((ETASQ+3.D0*RHSQ)*SN2HRT+
& (RHSQ+3.D0*ETASQ)*SN2ERT)
CUSVA(1,2) = 0.25D0*PNM* (H2* (RHSQ*SNERAT-ETASQ*SNHRAT)+
& 3.D0*SNSM)
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = 2.D0*FACL* ((RHSQ+3.D0*ETASQ)*SN2ERT-
& (ETASQ+3.D0*RHSQ)*SN2HRT)
END IF
C Compute A = U.Adjoint(V):
A(1,1) = U(1,1)*V(2,2) - U(1,2)*V(2,1)
A(1,2) = U(1,2)*V(1,1) - U(1,1)*V(1,2)
A(2,1) = A(1,2)
A(2,2) = V(1,1)*U(2,2) - U(2,1)*V(1,2)
C Compute T = (CU.A.Adjoint(CV) + SU.A.Adjoint(SV))*exp(-NUX-MUX):
IF (.NOT.SMALL) THEN
T(1,1) = A(1,1)*CNCM + 2.D0*A(1,2)* (NU*CMSN-MU*CNSM)/RAD +
& A(2,2)*SNSM/NM
T(1,2) = (A(1,1)*NM* (NU*CNSM-MU*CMSN)+
& A(1,2)* ((FAC* (NU2+MU2)+2.D0*NM*SNSM)* (NU2+MU2)-
& 4.D0*NM*NM*CNCM)/RAD+A(2,2)* (NU*CMSN-MU*CNSM))/RAD
T(2,1) = T(1,2)
T(2,2) = NM* (A(1,1)*SNSM+2.D0*A(1,2)* (NU*CNSM-MU*CMSN)/
& RAD) + A(2,2)*CNCM
ELSE
T(1,1) = A(1,1)*CNCM + 2.D0*P*A(1,2)*SUCVA(2,2) +
& A(2,2)*H2*SNHR1*SNHR2
T(1,2) = A(1,1)*PNM*NM*SUCVA(1,1) +
& 0.5D0*A(1,2)* (FAC+CNCM-H2*
& (ETASQ*SNHRAT+RHSQ*SNERAT)) + A(2,2)*P*SUCVA(2,2)
T(2,1) = T(1,2)
T(2,2) = NM* (A(1,1)*SNSM+2.D0*A(1,2)*PNM*SUCVA(1,1)) +
& A(2,2)*CNCM
END IF
C Compute U.Adjoint(V) = (CU.Uo+SU.Vo).Adjoint(SV.Uo+CV.Vo):
UVA(1,1) = DETU*CUSVA(1,1) + DETV*SUCVA(1,1) + T(1,1)
UVA(1,2) = DETU*CUSVA(1,2) + DETV*SUCVA(1,2) + T(1,2)
UVA(2,1) = DETU*CUSVA(2,1) + DETV*SUCVA(2,1) + T(2,1)
UVA(2,2) = DETU*CUSVA(2,2) + DETV*SUCVA(2,2) + T(2,2)
VUA(1,1) = UVA(2,2)
VUA(2,2) = UVA(1,1)
VUA(1,2) = -UVA(1,2)
VUA(2,1) = -UVA(2,1)
IF (.NOT.SMALL) THEN
DETCU = ((NU4+MU4)*CNCM-NM* ((NU2+MU2)*SNSM+2.D0*NM*FAC))/RAD2
DETCV = DETCU
DETSU = (2.D0*NM* (FAC-CNCM)+ (NU2+MU2)*SNSM)/ (PRAD2*NM)
DETSV = P*P*NM* (2.D0*NU3*MU3* (FAC-CNCM)+ (NU6+MU6)*SNSM)/
& RAD2
ELSE
DETCU = 0.25D0* (H2* (ETASQ*SNHRAT+RHSQ*SNERAT)+FAC+3.D0*CNCM)
DETCV = DETCU
RADH2 = ABS(RAD)*H2
AEHH = ABS(ETASQ-RHSQ)*H2
IF (RADH2.GE.0.01D0 .AND. AEHH.LT.0.01D0) THEN
FACL = H2*SNHR1*SNHR2
DETSU = (2.D0* (FAC-CNCM)+ (NU2+MU2)*FACL)/PRAD2
ELSE IF (RADH2.LT.0.01D0 .AND. AEHH.GE.0.01D0) THEN
FACL = H*0.5D0/P
DETSU = (SNERAT-SNHRAT)*FACL*FACL/ (ETASQ-RHSQ)
ELSE
FACL = 0.5D0*H2/P
DETSU = FACL*FACL*SNHDIF(AETAH,ARHH)*FAC
END IF
DETSV = (P*PNM/8.D0)* (2.D0*H*H*
& (RHSQ*RHSQ*SNERAT-ETASQ*ETASQ*SNHRAT)+
& 3.D0*NM* (CNCM-FAC)+15.D0* (ETASQ+RHSQ)*SNSM)
END IF
C Compute trace(CU.A.Adjoint(SU)).exp(-NUX-MUX):
BETA = A(1,1)*SUCVA(2,2) + 2.D0*A(1,2)*SUCVA(1,2) +
& A(2,2)*SUCVA(1,1)
C END IF
UUA(1,1) = DETU*DETCU + DETV*DETSU + BETA
UUA(1,2) = 0.D0
UUA(2,1) = 0.D0
UUA(2,2) = UUA(1,1)
C Compute trace(CV.A.Adjoint(SV)).exp(-NUX-MUX):
BETA = A(1,1)*CUSVA(2,2) + 2.D0*A(1,2)*CUSVA(1,2) +
& A(2,2)*CUSVA(1,1)
C END IF
VVA(1,1) = DETU*DETSV + DETV*DETCV + BETA
VVA(1,2) = 0.D0
VVA(2,1) = 0.D0
VVA(2,2) = VVA(1,1)
TRUVA = UVA(1,1) + UVA(2,2)
AC = UUA(1,1) - VVA(1,1)
B = -TRUVA
C = UUA(1,1) + VVA(1,1)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C Now we have got the value of M we can compute R.
ZCOSMT = COS(M/2.D0)
ZSINMT = SIN(M/2.D0)
DO 20 J = 1,2
DO 10 I = 1,2
R(I,J) = UT(I,J)*ZCOSMT - VT(I,J)*ZSINMT
10 CONTINUE
20 CONTINUE
C Compute the DET R and RINV
DETR = ZCOSMT**2*UUA(1,1) + ZSINMT**2*VVA(1,1) -
& ZSINMT*ZCOSMT* (UVA(1,1)+UVA(2,2))
C WRITE (8,FMT=*) 'DETR:',DETR
C
C WRITE (8,FMT=*) 'DETR*ADJUST-DET(R),NUX,MUX:',DETR*ADJUST-
C & R(1,1)*R(2,2)+R(1,2)*R(2,1),
C & NUX,MUX
C NUX = NUX-MUX
C NUX = 0.D0
NUX = MUX
C NUX = NUX + 2.D0*MUX
RINV(1,1) = R(2,2)/DETR
RINV(1,2) = -R(1,2)/DETR
RINV(2,1) = -R(2,1)/DETR
RINV(2,2) = R(1,1)/DETR
DO 40 J = 1,2
DO 30 I = 1,2
U(I,J) = (UUA(I,J)*ZCOSMT-UVA(I,J)*ZSINMT)/DETR
V(I,J) = (VUA(I,J)*ZCOSMT-VVA(I,J)*ZSINMT)/DETR
30 CONTINUE
40 CONTINUE
C That's done the very tricky renormalisation. Now do the SU and CV
C matrices which are the U and V matrices for Theta_{0}.
TRUVA = SUCVA(1,1) + SUCVA(2,2)
AC = (DETSU-DETCV)
B = -TRUVA
C = (DETSU+DETCV)
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C NOTE: in the following formulae we want CV and SU for BACKWARDS
C integration. We exploit the fact that the diagonal terms of
C SUCVA are odd functions of the direction of integration while
C the off-diagonal terms are even functions of the direction
C of integration, while the determinants which appear in the
C formulae are even functions of the direction of integration.
SU(1,1) = DETSU*COS(M/2.D0) + SUCVA(1,1)*SIN(M/2.D0)
SU(1,2) = -SUCVA(1,2)*SIN(M/2.D0)
SU(2,1) = -SUCVA(2,1)*SIN(M/2.D0)
SU(2,2) = DETSU*COS(M/2.D0) + SUCVA(2,2)*SIN(M/2.D0)
CV(1,1) = -SUCVA(2,2)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
CV(1,2) = -SUCVA(1,2)*COS(M/2.D0)
CV(2,1) = -SUCVA(2,1)*COS(M/2.D0)
CV(2,2) = -SUCVA(1,1)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
RETURN
END
C -------------------------------------------------------------------
SUBROUTINE NEGLM1(ELAM,NUX,RINV,P,S,Q,W,U,V,SU,CV,XO,XEND)
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,NUX,P,Q,S,W,XEND,XO
C ..
C .. Array Arguments ..
DOUBLE PRECISION CV(2,2),RINV(2,2),SU(2,2),U(2,2),V(2,2)
C ..
C .. Local Scalars ..
DOUBLE PRECISION AC,ALFA,ALFA2,ALFA4,AX,B,BETA,BETA2,BETA4,BX,C,
& CAXCBX,CAXSBX,COSBX,COSBX1,COSCOS,COSHAX,CSHCSH,
& D,DETCU,DETCV,DETR,DETSU,DETSV,DETU,DETV,FAC,H,M,
& PRAD,RAD,RAD2,RADA,RADB,RADM,SAXCBX,SAXSBX,SINBX,
& SINBX1,SINCOS,SINHAX,SINSIN,SNHCSH,SNHSNH,TR,
& TRUVA,ZCOSMT,ZSINMT
INTEGER I,J
C ..
C .. Local Arrays ..
DOUBLE PRECISION A(2,2),CU(2,2),CUSVA(2,2),R(2,2),SUCVA(2,2),
& SV(2,2),T(2,2),UT(2,2),UUA(2,2),UVA(2,2),VT(2,2),
& VUA(2,2),VVA(2,2)
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ATAN2,COS,EXP,SIGN,SIN,SQRT
C ..
C .. External Functions ..
DOUBLE PRECISION PHI
EXTERNAL PHI
C ..
C .. External Subroutines ..
EXTERNAL F06YAF
C ..
D = 2.D0*SQRT(P* (Q-ELAM*W))
ALFA = 0.5D0*SQRT((S+D)/P)
BETA = 0.5D0*SQRT((-S+D)/P)
DETU = U(1,1)*U(2,2) - U(2,1)*U(1,2)
DETV = V(1,1)*V(2,2) - V(2,1)*V(1,2)
H = XEND - XO
AX = ABS(ALFA*H)
FAC = EXP(-AX)
NUX = AX
IF (AX.GE.0.1D0) THEN
SINHAX = (1.D0-FAC*FAC)/2.D0
COSHAX = 1.D0 - SINHAX
SINHAX = SINHAX*SIGN(1.D0,H)/ALFA
ELSE
SINHAX = H*PHI(AX*AX)*FAC
COSHAX = (1.D0+FAC*FAC)/2.D0
END IF
BX = BETA*H
IF (ABS(BX).GE.0.1D0) THEN
SINBX1 = SIN(BX)/BETA
ELSE
SINBX1 = H*PHI(-BX*BX)
END IF
SINBX = SINBX1*FAC
COSBX1 = COS(BX)
COSBX = COSBX1*FAC
CAXCBX = COSHAX*COSBX1
CAXSBX = COSHAX*SINBX1
SAXSBX = SINHAX*SINBX1
SAXCBX = SINHAX*COSBX1
SNHCSH = (SINHAX*COSHAX)
SNHSNH = (SINHAX*SINHAX)
CSHCSH = (COSHAX*COSHAX)
SINCOS = (SINBX*COSBX)
SINSIN = (SINBX*SINBX)
COSCOS = (COSBX*COSBX)
ALFA2 = ALFA*ALFA
ALFA4 = ALFA2*ALFA2
BETA2 = BETA*BETA
BETA4 = BETA2*BETA2
RAD = ALFA2 + BETA2
RAD2 = RAD*RAD
PRAD = P*RAD
RADM = ALFA2 - BETA2
RADA = 3.D0*ALFA2 - BETA2
RADB = 3.D0*BETA2 - ALFA2
C Compute the CU,CV,SU,SV matrices here.
CU(1,1) = -RADM*SAXSBX*0.5D0 + CAXCBX
CU(1,2) = 0.5D0*SAXCBX + 0.5D0*CAXSBX
CU(2,1) = -0.5D0*RAD* (CAXSBX-SAXCBX)
CU(2,2) = CAXCBX + 0.5D0*RADM*SAXSBX
SU(1,1) = (SAXCBX-CAXSBX)/ (2.0D0*PRAD)
SU(1,2) = SAXSBX/ (2.0D0*P)
SU(2,1) = -SU(1,2)
SU(2,2) = (CAXSBX+SAXCBX)/ (2.0D0*P)
SV(1,1) = (((-2.D0*ALFA4+2.D0*BETA4+RAD2)*CAXSBX)+
& ((-2.D0*BETA4+2.D0*ALFA4+RAD2)*SAXCBX))* (0.5D0*P)
SV(1,2) = RAD2*SAXSBX*0.5D0*P
SV(2,1) = -SV(1,2)
SV(2,2) = (RADA*SAXCBX-RADB*CAXSBX)* (0.5D0*P)
CV(1,1) = -RADM*SAXSBX*0.5D0 + CAXCBX
CV(1,2) = (RAD*CAXSBX-RAD*SAXCBX)/2.0D0
CV(2,1) = -0.5D0* (CAXSBX+SAXCBX)
CV(2,2) = (RADM*SAXSBX*0.5D0) + CAXCBX
C Now we can compute UT and VT
CALL F06YAF('N','N',2,2,2,1.D0,CU,2,U,2,0.D0,UT,2)
CALL F06YAF('N','N',2,2,2,1.D0,SU,2,V,2,1.D0,UT,2)
CALL F06YAF('N','N',2,2,2,1.D0,SV,2,U,2,0.D0,VT,2)
CALL F06YAF('N','N',2,2,2,1.D0,CV,2,V,2,1.D0,VT,2)
C Compute SU.Adjoint(CV).exp(-AX):
SUCVA(1,1) = 0.5D0* (SNHCSH-SINCOS)/PRAD
SUCVA(1,2) = 0.25D0* (SNHSNH+SINSIN)/P
SUCVA(2,1) = SUCVA(1,2)
SUCVA(2,2) = 0.5D0* (SNHCSH+SINCOS)/P
C Compute CU.Adjoint(SV).exp(-AX):
CUSVA(1,1) = 0.5D0* (RADA*SNHCSH-RADB*SINCOS)/P
CUSVA(1,2) = -0.25D0*RAD* (SNHSNH*BETA2-3.D0*CSHCSH+SINSIN*ALFA2+
& 3.D0*COSCOS)/P
CUSVA(2,1) = CUSVA(1,2)
CUSVA(2,2) = 0.5D0*PRAD* (RADA*SNHCSH+RADB*SINCOS)
C Compute A = U.Adjoint(V)(xo):
A(1,1) = U(1,1)*V(2,2) - U(1,2)*V(2,1)
A(1,2) = U(1,2)*V(1,1) - U(1,1)*V(1,2)
A(2,1) = A(1,2)
A(2,2) = U(2,2)*V(1,1) - U(2,1)*V(1,2)
C Compute T = (CU.A.Adjoint(CV) + SU.Adjoint(A).Adjoint(SV))*exp(-AX):
T(1,1) = A(1,1)* (CSHCSH-SINSIN*BETA2) + A(1,2)* (SNHCSH+SINCOS) +
& A(2,2)* (CSHCSH-COSCOS)/RAD
T(1,2) = 0.5D0*RAD*A(1,1)* (SNHCSH-SINCOS) +
& 0.5D0*A(1,2)* (CSHCSH+COSCOS+BETA2*SNHSNH-ALFA2*SINSIN) +
& 0.5D0*A(2,2)* (SNHCSH+SINCOS)
T(2,1) = T(1,2)
T(2,2) = A(1,1)*RAD* (CSHCSH-COSCOS) +
& A(1,2)*RAD* (SNHCSH-SINCOS) +
& A(2,2)* (CSHCSH-SINSIN*BETA2)
C Compute U.Adjoint(V)(xend) = (CU.Uo+SU.Vo).Adjoint(SV.Uo+CV.Vo):
UVA(1,1) = DETU*CUSVA(1,1) + DETV*SUCVA(1,1) + T(1,1)
UVA(1,2) = DETU*CUSVA(1,2) + DETV*SUCVA(1,2) + T(1,2)
UVA(2,1) = DETU*CUSVA(2,1) + DETV*SUCVA(2,1) + T(2,1)
UVA(2,2) = DETU*CUSVA(2,2) + DETV*SUCVA(2,2) + T(2,2)
VUA(1,1) = UVA(2,2)
VUA(2,2) = UVA(1,1)
VUA(1,2) = -UVA(1,2)
VUA(2,1) = -UVA(2,1)
DETCU = 0.25D0*RADM* (SNHSNH+SINSIN) + 0.5D0* (CSHCSH+COSCOS)
DETSU = 0.25* (SNHSNH-SINSIN)/ (P*PRAD)
DETCV = DETCU
DETSV = 0.25D0*P*PRAD* (RADA*RADA*SNHSNH-RADB*RADB*SINSIN)
C Compute trace(CU.A.Adjoint(SU)).exp(-AX):
TR = 0.5D0* (A(1,1)* (SNHCSH+SINCOS)+A(1,2)* (SNHSNH+SINSIN))/P +
& 0.5D0*A(2,2)* (SNHCSH-SINCOS)/PRAD
UUA(1,1) = DETU*DETCU + DETV*DETSU + TR
UUA(1,2) = 0.D0
UUA(2,1) = 0.D0
UUA(2,2) = UUA(1,1)
C Compute trace(CV.Adjoint(A).Adjoint(SU)).exp(-AX):
TR = 0.5D0*P*RADA* (RAD* (A(1,1)*SNHCSH+A(1,2)*SNHSNH)+
& A(2,2)*SNHCSH) + 0.5D0*P*RADB*
& (RAD* (A(1,1)*SINCOS+A(1,2)*SINSIN)-A(2,2)*SINCOS)
VVA(1,1) = DETU*DETSV + DETV*DETCV + TR
VVA(1,2) = 0.D0
VVA(2,1) = 0.D0
VVA(2,2) = VVA(1,1)
TRUVA = UVA(1,1) + UVA(2,2)
AC = UUA(1,1) - VVA(1,1)
B = -TRUVA
C = UUA(1,1) + VVA(1,1)
C We now choose M to maximise the function
C
C ABS(AC.cos(M) + BC.sin(M) + C):
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
C Now we have got the value of M we can compute R.
ZCOSMT = COS(M/2.D0)
ZSINMT = SIN(M/2.D0)
DO 20 J = 1,2
DO 10 I = 1,2
R(I,J) = (UT(I,J)*ZCOSMT-VT(I,J)*ZSINMT)
10 CONTINUE
20 CONTINUE
C Compute the DET R and RINV
DETR = ZCOSMT**2*UUA(1,1) + ZSINMT**2*VVA(1,1) -
& ZSINMT*ZCOSMT* (UVA(1,1)+UVA(2,2))
RINV(1,1) = R(2,2)/DETR
RINV(1,2) = -R(1,2)/DETR
RINV(2,1) = -R(2,1)/DETR
RINV(2,2) = R(1,1)/DETR
U(1,1) = (UUA(1,1)*ZCOSMT-UVA(1,1)*ZSINMT)/DETR
U(1,2) = (-UVA(1,2)*ZSINMT)/DETR
U(2,1) = (-UVA(2,1)*ZSINMT)/DETR
U(2,2) = (UUA(2,2)*ZCOSMT-UVA(2,2)*ZSINMT)/DETR
V(1,1) = (VUA(1,1)*ZCOSMT-VVA(1,1)*ZSINMT)/DETR
V(1,2) = (VUA(1,2)*ZCOSMT)/DETR
V(2,1) = (VUA(2,1)*ZCOSMT)/DETR
V(2,2) = (VUA(2,2)*ZCOSMT-VVA(2,2)*ZSINMT)/DETR
C That's done the very tricky renormalisation. Now do the SU and CV
C matrices which are the U and V matrices for Theta_{0}.
TRUVA = SUCVA(1,1) + SUCVA(2,2)
AC = (DETSU-DETCV)
B = -TRUVA
C = (DETSU+DETCV)
IF (C.GE.0.D0) THEN
M = ATAN2(B,AC)
ELSE
M = ATAN2(-B,-AC)
END IF
SU(1,1) = DETSU*COS(M/2.D0) + SUCVA(1,1)*SIN(M/2.D0)
SU(1,2) = -SUCVA(1,2)*SIN(M/2.D0)
SU(2,1) = -SUCVA(2,1)*SIN(M/2.D0)
SU(2,2) = DETSU*COS(M/2.D0) + SUCVA(2,2)*SIN(M/2.D0)
CV(1,1) = -SUCVA(2,2)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
CV(1,2) = -SUCVA(1,2)*COS(M/2.D0)
CV(2,1) = -SUCVA(2,1)*COS(M/2.D0)
CV(2,2) = -SUCVA(1,1)*COS(M/2.D0) - DETCV*SIN(M/2.D0)
RETURN
END
C =====================================================================
C ==================== SUBROUTINES FOR AUTOMATIC MESHING ==============
C =====================================================================
SUBROUTINE MESH4(XO,XEND,ELAM,SL4COF,XMESH,NMESH,IC,IMESH,TOL,
& YSTOR1,YSTOR2,MULTI,IM1ST,CNMESH,WORK,IWORK,
& IFAIL)
C
C The purpose of this subroutine is to compute a mesh which can be
C used as an initial mesh for solving a 4th-order Sturm-Liouville
C eigenvalue problems.
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,TOL,XEND,XO
INTEGER CNMESH,IC,IFAIL,IM1ST,IMESH,IWORK,MULTI,NMESH
C ..
C .. Array Arguments ..
DOUBLE PRECISION WORK(0:IWORK,1:5),XMESH(0:NMESH),
& YSTOR1(4,0:CNMESH),YSTOR2(4,0:CNMESH)
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4COF
C ..
C .. Parameters ..
DOUBLE PRECISION HLF,ONE,SAFE,TRD,FIVE,SAFE3,P008
PARAMETER (HLF=5.d-1,ONE=1.d0,SAFE=0.8D0,TRD=ONE/3.d0,FIVE=5.d0,
& SAFE3=7.29d-1,P008=8.d-3)
C ..
C .. Local Scalars ..
DOUBLE PRECISION DI,EPS,ERREST,H,HMAX,RATIO,TEND,TO,TOLOC,TOLOW,
& X1,X3
INTEGER I,ICNT,IDI,IRF,J,KNTR
LOGICAL PHASE1
C ..
C .. External Functions ..
DOUBLE PRECISION GETERR,X02AJF
EXTERNAL GETERR,X02AJF
C ..
C .. External Subroutines ..
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,INT,LOG,MAX,MIN,SIGN
C ..
TO = XO/ (1.D0+ABS(XO))
TEND = XEND/ (1.D0+ABS(XEND))
IFAIL = 0
EPS = X02AJF(0.D0)*1.D2
IF (TO.EQ.TEND) THEN
IMESH = 0
XMESH(IC) = XO
RETURN
ELSE
DI = SIGN(1.D0,TEND-TO)
IDI = DI
END IF
IF (IDI.GT.0) THEN
I = IC
ELSE
I = NMESH
END IF
XMESH(I) = TO
ICNT = 0
KNTR = 0
IRF = INT(LOG(X02AJF(0.D0))*HLF/LOG(SAFE))
TOLOC = TOL
TOLOW = TOLOC*SAFE
PHASE1 = .TRUE.
C Set initial stepsize:
H = MIN(ABS(TEND-TO),1.D0)/20.D0
C Set maximum stepsize:
HMAX = ABS(TEND-TO)/20.D0
10 X3 = XMESH(I) + H*DI
X1 = XMESH(I)
IF (IDI.GT.0) X3 = MIN(X3,TEND)
IF (IDI.LT.0) X3 = MAX(X3,TEND)
C Get an error estimate.
ERREST = GETERR(X1,X3,SL4COF,ELAM,YSTOR1,YSTOR2,MULTI,IM1ST,
& CNMESH,WORK,IWORK)
IF (ERREST.GT.TOLOC) THEN
RATIO = SAFE
IF (TOLOC.LT.ERREST*SAFE3) RATIO = (TOLOC/ERREST)**TRD
H = H*RATIO
ICNT = ICNT + 1
IF (ICNT.LT.IRF) GO TO 10
IFAIL = 10
RETURN
END IF
IF (ERREST.LT.TOLOW .AND. ICNT.EQ.0) THEN
RATIO = FIVE
IF (ERREST.GT.TOLOW*P008) RATIO = (TOLOW/ERREST)**TRD
H = MIN(HMAX,H*RATIO)
IF (PHASE1 .AND. H.LT.HMAX) THEN
KNTR = KNTR + 1
IF (KNTR.LT.5) GO TO 10
END IF
END IF
PHASE1 = .FALSE.
ICNT = 0
I = I + IDI
XMESH(I) = X3
IF ((IDI.GT.0.AND.X3.LT.TEND-EPS) .OR.
& (IDI.LT.0.AND.X3.GT.TEND+EPS)) THEN
C We have not yet reached the end of the range. Meshing must continue
C so we must check that there is space to hold at least another one
C mesh-point.
IF (IDI.GT.0 .AND. I.GE.NMESH) THEN
IFAIL = 11
RETURN
END IF
IF (IDI.LT.0 .AND. I.LE.IC) THEN
IFAIL = 11
RETURN
END IF
GO TO 10
END IF
C If we are here, it means that the meshing has finished. We must
C move the mesh into the place where the calling routine expects
C to find it.
IF (IDI.LT.0 .AND. I.GT.IC) THEN
DO 20 J = 0,NMESH - I
XMESH(IC+J) = XMESH(I+J)
20 CONTINUE
END IF
IF (IDI.GT.0) THEN
IMESH = I - IC
C Transform mesh back to original coordinates
DO 30 J = IC,I
XMESH(J) = XMESH(J)/ (1.D0-ABS(XMESH(J)))
30 CONTINUE
ELSE
IMESH = NMESH - I
C Transform mesh back to original coordinates
DO 40 J = IC,IMESH + IC
XMESH(J) = XMESH(J)/ (1.D0-ABS(XMESH(J)))
40 CONTINUE
END IF
RETURN
END
C ---------------------------------------------------------------------
DOUBLE PRECISION FUNCTION GETERR(X1,X3,SL4COF,ELAM,YSTOR1,YSTOR2,
& MULTI,IMATCH,NMESH,WORK,IWORK)
C .. External Subroutines ..
EXTERNAL EFN
C ..
C .. Local Arrays ..
DOUBLE PRECISION Y(1:4)
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION ELAM,X1,X3
INTEGER IMATCH,IWORK,MULTI,NMESH
C ..
C .. Array Arguments ..
DOUBLE PRECISION WORK(0:IWORK,5),YSTOR1(4,0:NMESH),
& YSTOR2(4,0:NMESH)
C ..
C .. Local Scalars ..
DOUBLE PRECISION DDPY21,DDPY23,DY21,DY23,ERRST,FAC1,FAC3,P1,P2,P3,
& Q1,Q2,Q3,S1,S2,S3,W1,W2,W3,X2,Y21,Y23
LOGICAL HOT
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,MAX
C ..
C .. Subroutine Arguments ..
EXTERNAL SL4COF
C ..
X2 = (X3+X1)*0.5D0
C At the moment we only form the error estimate using one of the
C eigenfunctions, in the case of an eigenvalue of multiplicity 2.
FAC1 = 1.D0/ (1.D0-ABS(X1))**2
FAC3 = 1.D0/ (1.D0-ABS(X3))**2
CALL SL4COF(X1/ (1.D0-ABS(X1)),P1,S1,Q1,W1)
CALL SL4COF(X2/ (1.D0-ABS(X2)),P2,S2,Q2,W2)
CALL SL4COF(X3/ (1.D0-ABS(X3)),P3,S3,Q3,W3)
HOT = .FALSE.
CALL EFN(WORK(0,1),WORK(1,2),WORK(1,3),WORK(1,4),WORK(1,5),ELAM,
& X1/ (1.D0-ABS(X1)),YSTOR1,IMATCH,NMESH,Y,HOT)
Y21 = Y(1)**2
DY21 = Y(2)**2
DDPY21 = Y(4)**2
Y21 = MAX(Y21,ABS(Y(1)))
DY21 = MAX(DY21,ABS(Y(2)))
DDPY21 = MAX(DDPY21,ABS(Y(4)))
CALL EFN(WORK(0,1),WORK(1,2),WORK(1,3),WORK(1,4),WORK(1,5),ELAM,
& X3/ (1.D0-ABS(X3)),YSTOR1,IMATCH,NMESH,Y,HOT)
Y23 = Y(1)**2
DY23 = Y(2)**2
DDPY23 = Y(4)**2
Y23 = MAX(Y23,ABS(Y(1)))
DY23 = MAX(DY23,ABS(Y(2)))
DDPY23 = MAX(DDPY23,ABS(Y(4)))
ERRST = ((Q1-Q2)-ELAM* (W1-W2))*Y21 + (S1-S2)*DY21 +
& ((1.D0/P2)- (1.D0/P1))*DDPY21
ERRST = FAC3* (((Q3-Q2)-ELAM* (W3-W2))*Y23+ (S3-S2)*DY23+
& ((1.D0/P2)- (1.D0/P3))*DDPY23) + FAC1*ERRST
ERRST = ERRST* ((X3-X1)/6.D0)
GETERR = ABS(ERRST)
RETURN
END
C ---------------------------------------------------------------
SUBROUTINE COARS4(N,A,B,XMESH,ISING)
C .. Scalar Arguments ..
DOUBLE PRECISION A,B
INTEGER ISING,N
C ..
C .. Array Arguments ..
DOUBLE PRECISION XMESH(0:N)
C ..
C .. Local Scalars ..
DOUBLE PRECISION H
INTEGER I
C ..
C .. Intrinsic Functions ..
INTRINSIC DBLE
C ..
XMESH(0) = A
XMESH(N) = B
IF (ISING.EQ.0) THEN
C THE PROBLEM IS REGULAR, USE A UNIFORM MESH
H = (B-A)/DBLE(N)
DO 10 I = 1,N - 1
XMESH(I) = XMESH(I-1) + H
10 CONTINUE
ELSE IF (ISING.EQ.1) THEN
C X = A IS SINGULAR BUT X = B IS NOT
H = (B-A)/DBLE(N-4)
XMESH(1) = XMESH(0) + H/16.D0
XMESH(2) = XMESH(1) + H/16.D0
XMESH(3) = XMESH(2) + H/8.D0
XMESH(4) = XMESH(3) + H/4.D0
XMESH(5) = XMESH(4) + H/2.D0
DO 20 I = 6,N - 1
XMESH(I) = XMESH(I-1) + H
20 CONTINUE
ELSE IF (ISING.EQ.2) THEN
C X=A IS REGULAR BUT X = B IS SINGULAR
H = (B-A)/DBLE(N-4)
XMESH(N-1) = XMESH(N) - H/16.D0
XMESH(N-2) = XMESH(N-1) - H/16.D0
XMESH(N-3) = XMESH(N-2) - H/8.D0
XMESH(N-4) = XMESH(N-3) - H/4.D0
XMESH(N-5) = XMESH(N-4) - H/2.D0
DO 30 I = N - 6,1,-1
XMESH(I) = XMESH(I+1) - H
30 CONTINUE
ELSE IF (ISING.EQ.3) THEN
C BOTH ENDS OF THE INTERVAL ARE SINGULAR POINTS.
H = (B-A)/DBLE(N-8)
XMESH(1) = XMESH(0) + H/16.D0
XMESH(2) = XMESH(1) + H/16.D0
XMESH(3) = XMESH(2) + H/8.D0
XMESH(4) = XMESH(3) + H/4.D0
XMESH(5) = XMESH(4) + H/2.D0
XMESH(N-1) = XMESH(N) - H/16.D0
XMESH(N-2) = XMESH(N-1) - H/16.D0
XMESH(N-3) = XMESH(N-2) - H/8.D0
XMESH(N-4) = XMESH(N-3) - H/4.D0
XMESH(N-5) = XMESH(N-4) - H/2.D0
DO 40 I = 6,N - 6
XMESH(I) = XMESH(I-1) + H
40 CONTINUE
ELSE
C THERE IS AN ERROR IN THE VALUE OF ISING
END IF
RETURN
END
C =====================================================================
C ========================== AUXILIARY SUBROUTINES ====================
C =====================================================================
DOUBLE PRECISION FUNCTION PHI(V)
C .. Parameters ..
DOUBLE PRECISION A0,A1,A2,A3,A4,A5
PARAMETER (A0=1.D0,A1=A0/6.D0,A2=A1/20.D0,A3=A2/42.D0,A4=A3/72.D0,
& A5=A4/110.D0)
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION V
C ..
PHI = A0 + V* (A1+V* (A2+V* (A3+V* (A4+V*A5))))
RETURN
END
C ---------------------------------------------------------------------
DOUBLE PRECISION FUNCTION SNHDIF(ETA,H)
C This function computes the quantity
C (sinh(eta)/eta - sinh(h)/h)/(eta**2-h**2)
C in the case where eta and h are both small.
C
C .. Parameters ..
DOUBLE PRECISION FAC1,FAC2,FAC3,FAC4
PARAMETER (FAC1=1.D0/6.D0,FAC2=FAC1/20.D0,FAC3=FAC2/42.D0,
& FAC4=FAC3/72.D0)
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION ETA,H
C ..
C .. Local Scalars ..
DOUBLE PRECISION ETA2,ETA4,ETA6,H2,H4,H6
C ..
H2 = H*H
H4 = H2*H2
H6 = H4*H2
ETA2 = ETA*ETA
ETA4 = ETA2*ETA2
ETA6 = ETA4*ETA2
SNHDIF = FAC1 + FAC2* (ETA2+H2) + FAC3* (ETA4+H4+ETA2*H2) +
& FAC4* (ETA6+H6+ETA4*H2+ETA2*H4)
RETURN
END
C ---------------------------------------------------------------------
SUBROUTINE MONIT(A,IA,N,IPRN)
C .. Scalar Arguments ..
INTEGER IA,IPRN,N
C ..
C .. Array Arguments ..
DOUBLE PRECISION A(IA,N)
C ..
C .. Local Scalars ..
INTEGER I,J
C ..
C .. Intrinsic Functions ..
INTRINSIC REAL
C ..
DO 10 J = 1,N
WRITE (IPRN,FMT=*) (REAL(A(J,I)),I=1,N)
10 CONTINUE
RETURN
END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
# End of shell archive
exit 0
.