PROGRAM RXP DEMO
C This program demonstrates the use of the coefficients calculated
C by the programs axitrx. The program will produce a table of values
C for different values of lambda and s as set in the DATA statements.
INTEGER MAXS,MAXL
PARAMETER(MAXS=5,MAXL=6)
DOUBLE PRECISION DLAMDA(MAXL),S(MAXS),RXP,VAL(MAXL)
INTEGER I,J
DATA DLAMDA/1D0,2D0,3D0,4D0,5D0,0.5D0/
DATA S/2.01D0,2.02D0,2.05D0,2.1D0,2.5D0/
WRITE(6,10)
10 FORMAT(1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XP11'/)
WRITE(6,20)
20 FORMAT(' s',6X,25('-'),' lambda ',24('-'),/, 14X,
& '1 2 3 4 5 1/2')
DO 50 I=1,MAXS
DO 40 J=1,MAXL
40 VAL(J)=RXP(1,1,S(I),DLAMDA(J))
50 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
WRITE(6,110)
110 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XP12'/)
WRITE(6,20)
DO 150 I=1,MAXS
DO 140 J=1,MAXL
140 VAL(J)=RXP(1,2,S(I),DLAMDA(J))
150 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
WRITE(6,210)
210 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XP21'/)
WRITE(6,20)
DO 250 I=1,MAXS
DO 240 J=1,MAXL
240 VAL(J)=RXP(2,1,S(I),DLAMDA(J))
250 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
WRITE(6,310)
310 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XP22'/)
WRITE(6,20)
DO 350 I=1,MAXS
DO 340 J=1,MAXL
340 VAL(J)=RXP(2,2,S(I),DLAMDA(J))
350 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
1100 FORMAT(F8.3,8F10.4)
STOP
END
DOUBLE PRECISION FUNCTION RXP(IALPHA,IBETA,S,DLAMDA)
C -------------------- Arguments of function ---------------------------
DOUBLE PRECISION DLAMDA
INTEGER IALPHA,IBETA
C ----------------------------------------------------------------------
C P
C This routine calculates the functions X (s,lambda)
C ALPHA BETA
C defined in Jeffrey, Phys. Fluids(1992).
C Coefficients calculated by one of the programs AXITRx or AXISTx
C are summed using equations (20a,20b).
C The coefficients are expected in RXP200.DAT in the format on the
C distribution disk.
C -------------------------- Local variables -----------------------------
DOUBLE PRECISION RXPCOF(0:45450)
DOUBLE PRECISION G1,G2,G3,G1SUM,G2SUM,G3SUM
DOUBLE PRECISION T,ZT,ZL,COEF,EPS,S
DOUBLE PRECISION XL,XLP1,XL13, XL12
INTEGER I,J,ISTEP,IQ,MMIN,MMAX,M,NM,ISIGN
LOGICAL LOADED
CHARACTER*4 ID
DATA LOADED/.FALSE./
C ----------------------- Start of executable code ---------------------
IF(.NOT.LOADED) THEN
OPEN(1,FILE='rxp300.dat',STATUS='OLD')
READ(1,10) ID,MAXS
10 FORMAT(A4,I5)
IF(ID.NE.' RXP') THEN
WRITE(*,20)
20 FORMAT(1X,'WRONG COEFFICIENTS FOR RYC')
STOP
ENDIF
READ(1,30) RXPCOF
30 FORMAT(D22.16)
CLOSE(1)
LOADED=.TRUE.
ENDIF
IF(DLAMDA.GT.1D0) THEN
I=3-IALPHA
J=3-IBETA
XL=1D0/DLAMDA
C RXP and RGX change sign when we invert lambda
ISIGN = -1
ELSE
I=IALPHA
J=IBETA
XL=DLAMDA
ISIGN = 1
ENDIF
XLP1=XL+1D0
XL12=XLP1**2
XL13=XLP1**3
C These are the g_i for X^G_{11,12} (i=1..3 in paper)
IF(I.EQ.1) THEN
G1=3D0*XL**2/XL13
G2=3d0*(XL-4D0*XL**2)/(10D0*XL12)
G3=(5D0-97D0*XL+64D0*XL**2-44D0*XL**3+XL**4)/(140D0*XL12)
ISTEP=1
ELSE
C These are the g_i for X^G_{21,22}
G1=3D0*XL/XL13
G2=3d0*(XL-4D0)/(10D0*XL12)
G3=(5D0*XL**4-97D0*XL**3+64D0*XL**2-44D0*XL+1D0)
& /(140D0*XL12*XL**2)
C These expressions essentially invert lambda.
ISIGN=-ISIGN
ISTEP=-1
ENDIF
T=2D0/S
EPS=0.25D0*S*S-1D0
IF(I.NE.J) THEN
ZT=T*T
C Notice that because the first coef is a signature we supply the 0.
RXP=0D0+G1/EPS-(G2+G3*EPS)*DLOG( 1D0-4D0/(S*S) )-G3
MMIN=2
MMAX=MAXS
ELSE
ZT=T
RXP=G1*S/(2D0*EPS)+(G2+G3*EPS)*DLOG((S+2D0)/(S-2D0))-G3*S
MMIN=1
MMAX=MAXS-1
ENDIF
T=T*T
G1SUM=G1
G2SUM=2D0*G2
G3SUM=4D0*G3
DO 100 M=MMIN,MMAX,2
IF (ISTEP.GT.0) THEN
NM = (M*(M+1))/2
ELSE
NM = (M*(M+1))/2 +M+2
ENDIF
ZL=XLP1**(-M)
COEF=0D0
C Some computers will underflow in the next loop.
DO 50 IQ=0,M
NM=NM+ISTEP
COEF=COEF+RXPCOF(NM)*ZL
ZL=ZL*XL
50 CONTINUE
RXP=RXP+ZT*(COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) )
ZT=ZT*T
100 CONTINUE
RXP = DBLE(ISIGN)*RXP
IF(I.EQ.J)RETURN
RXP=-RXP*4D0/XLP1**2
IF(I.EQ.2)RXP=RXP*XL**2
RETURN
END
C
BLOCK DATA RXPINI
DOUBLE PRECISION RXPCOF(45450)
COMMON /XPD250/RXPCOF
DATA RXPCOF(1)/0D0/
END