PROGRAM RXQ DEM C This program demonstrates the use of the coefficients calculated C by the programs axistx. The program will produce a table of values C for different values of LAMBDA and S as set in the DATA statements. C INTEGER MAXS,MAXL PARAMETER(MAXS=5,MAXL=6) DOUBLE PRECISION DLAMDA(MAXL),S(MAXS),RXQ,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 XQ11'/) 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)=RXQ(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 XQ12'/) WRITE(6,20) DO 150 I=1,MAXS DO 140 J=1,MAXL 140 VAL(J)=RXQ(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 XQ21'/) WRITE(6,20) DO 250 I=1,MAXS DO 240 J=1,MAXL 240 VAL(J)=RXQ(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 XQ22'/) WRITE(6,20) DO 350 I=1,MAXS DO 340 J=1,MAXL 340 VAL(J)=RXQ(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 C DOUBLE PRECISION FUNCTION RXQ(IALPHA,IBETA,S,DLAMDA) C ---------------------- Arguments of function ------------------------ INTEGER IALPHA, IBETA DOUBLE PRECISION S, DLAMDA C --------------------------------------------------------------------- C Q C This routine calculates the functions X (s,lambda) C ALPHA BETA C defined in Jeffrey, ' C ',Phys. Fluids (1991). C Coefficients calculated by one of the programs AXISTx are summed. C It goes forwards or backwards through the coefficients ( ISTEP) C depending upon whether I=1 or 2. C This technique is not explained in the paper. C Starting at equn (5.d) we see that X can be obtained from X C 22 11 C by inverting lambda. Now looking at equation (47) we substitute C 1/lambda and redefine q as k-q. Factoring the lambda**k, we see that C lambda**q now multiplies P(1,q,k-q), which means running 'backwards' C through the coefficients relative to the published case. C C The program reads the coefficients from RXQ200.DAT only once and checks C to make certain the signature is correct (cf AXIST1). C C ------------------------ Local variables ---------------------------- DOUBLE PRECISION RXQCOF(0:45450) DOUBLE PRECISION G1, G2, G3, G1SUM, G2SUM, G3SUM DOUBLE PRECISION T, ZT, ZL, COEF, EPS, XL, XLP1, XL13 INTEGER I, J, ISTEP, IQ, MMIN, MMAX, M, NM LOGICAL LOADED CHARACTER*4 ID DATA LOADED/.FALSE./ C ------------------------ Start of executable code ------------------- IF(.NOT.LOADED) THEN OPEN(1,FILE='rxq300.dat',STATUS='OLD') READ(1,10) ID,MAXS 10 FORMAT(A4,I5) IF(ID.NE.' RXQ') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RXQ') STOP ENDIF READ(1,30) RXQCOF 30 FORMAT(D22.16) CLOSE(1) LOADED=.TRUE. ENDIF IF(DLAMDA.GT.1D0) THEN I=3-IALPHA J=3-IBETA XL=1D0/DLAMDA ELSE I=IALPHA J=IBETA XL=DLAMDA ENDIF XLP1=XL+1D0 XL13=XLP1**3 IF(I.EQ.J) THEN IF(I.EQ.1) THEN C These are the g_i for X^Q_{11} (i=1..3 in paper) C this means lambda <1 and i=k=1. ie xq11, ie even G1=3D0*XL**2/(2D0*XL13) G2=3D0*(XL+2D0*XL**2-9D0*XL**3)/(20D0*XL13) G3=(5D0-XL-201D0*XL**2+251D0*XL**3-184D0*XL**4) & /(280D0*XL13) ISTEP=1 ELSE C These are the g_i for X^Q_{22} G1=3D0*XL/(2D0*XL13) G2=3D0*(XL**2+2D0*XL-9D0)/(20D0*XL13) G3=(5D0*XL**4-XL**3-201D0*XL**2+251D0*XL-184D0) & /(280D0*XL*XL13) ISTEP=-1 ENDIF ELSE IF(I.EQ.1) THEN C These are the g_i for X^Q_{12} (i=4..6 in paper) G1=3D0*XL**3/(2D0*XL13) G2=3D0*(-2D0*XL**2+1D0*XL**3-2D0*XL**4)/(10D0*XL13) G3=(-65D0*XL+34D0*XL**2-411D0*XL**3+76D0*XL**4 & -44D0*XL**5)/(280D0*XL13) ISTEP=1 ELSE C These are the g_i for X^M_21 G1=3D0/(2D0*XL13) G2=3D0*(-2D0*XL**2+1D0*XL-2D0)/(10D0*xl*XL13) G3=(-65D0*XL**4+34D0*XL**3-411D0*XL**2+76D0*XL & -44D0)/(280D0*XL13*xl**2) ISTEP=-1 ENDIF ENDIF T=2D0/S EPS=0.25D0*S*S-1D0 IF(I.EQ.J) THEN ZT=T*T C Notice that because the first coef is a signature we supply the 1. RXQ=1D0+G1/EPS+(G2+G3*EPS)*DLOG( (S*S)/(S*S-4D0))-G3 MMIN=2 MMAX=200 ELSE ZT=T RXQ=G1*S/(2D0*EPS)+(G2+G3*EPS)*DLOG((S+2D0)/(S-2D0))-G3*S MMIN=1 MMAX=199 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 ENDIF C This loop adds up the powers of lambda to obtain the term f_m in the C paper. Notice that the 2**(-m) is not needed, because the coefficients C in the file already contain that factor. ZL=1D0 COEF=0D0 C Some computers will underflow in the next loop. DO 50 IQ=0,M COEF=COEF+RXQCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE COEF=COEF*XLP1**(-M) RXQ=RXQ + (COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) )*ZT ZT=ZT*T 100 CONTINUE IF(I.EQ.J)RETURN RXQ=RXQ*8D0/XLP1**3 IF(I.EQ.2)RXQ=RXQ*XL**3 RETURN END C BLOCK DATA RXQINI DOUBLE PRECISION RXQCOF(0:45450) COMMON /XQD200/RXQCOF DATA RXQCOF(1)/0D0/ END