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