PROGRAM RYB DEM C This program demonstrates the use of the coefficients calculated C by either the programs trntrx or trnrox. The program will produce C a table of values for different values of lambda and s as set in C the DATA statements. C DOUBLE PRECISION DLAMDA(5),S(4),RYB,RYB11,RYB12,RYB21,RYB22 INTEGER I,J DATA DLAMDA/1D0,.5D0,.25D0,.125D0,2D0/ DATA S/2.01D0,2.1D0,2.5D0,3D0/ WRITE(6,200) 200 FORMAT(1X,'Table of values for Resistance Functions YB') DO 50 J=1,5 WRITE(6,1000)DLAMDA(J) 1000 FORMAT(/1X,'LAMBDA =',F8.3) WRITE(6,100) 100 FORMAT(/10X,'S',10X,'RYB(1,1)',8X,'RYB(1,2)', & 8X,'RYB(2,1)',8X,'RYB(2,2)') DO 50 I=1,4 RYB11=RYB(1,1,S(I),DLAMDA(J)) RYB12=RYB(1,2,S(I),DLAMDA(J)) RYB21=RYB(2,1,S(I),DLAMDA(J)) RYB22=RYB(2,2,S(I),DLAMDA(J)) 50 WRITE(6,1100)S(I),RYB11,RYB12,RYB21,RYB22 1100 FORMAT(F13.3,4F16.5) STOP END C DOUBLE PRECISION FUNCTION RYB(IALPHA,IBETA,S,DLAMDA) C ---------------------- Arguments of function -------------------------- DOUBLE PRECISION S,DLAMDA INTEGER IALPHA,IBETA C ----------------------------------------------------------------------- C B C This routine calculates the functions Y (s,lambda) C ALPHA BETA C defined in Jeffrey & Onishi (1984) J. Fluid Mech. Vol 139, 261-290. C Coefficients calculated by one of the programs TRNTRx are summed using C equations (5.9)-(5.10). The subprogram uses equation (1.9a) to convert C any call to one in which lambda.LT.1. It goes forwards or backwards C through the coefficients ( ISTEP) depending upon whether I=1 or 2. C This technique is not explained in the paper. C Starting at equn 1.9 we see that Y can be obtained from Y C 22 11 C by inverting lambda. Now looking at equn 3.15 we substitute 1/lambda and C redefine q as k-q. Factoring the lambda**k, we see that lambda**q now C multiplies P(1,q,k-q), which means running 'backwards' through the C coefficients relative to the published case. C C In the program the coefficients are read in and checked C to make certain the signature is correct (cf TRNTR1). C The coefficients are expected in RYBxxx.DAT C C ------------------------ Local variables ----------------------------- DOUBLE PRECISION RYBCOF(0:45450),XL,XLP1,G2,G3,G2SUM,G3SUM DOUBLE PRECISION T,ZT,ZL,COEF,EPS INTEGER I,J,ISIGN,ISTEP,IQ,MMIN,MMAX,M,NM, MAXS LOGICAL LOADED CHARACTER*4 ID DATA LOADED/.FALSE./ C ----------------------- Start of executable code --------------------- IF(.NOT.LOADED) THEN OPEN(3,FILE='ryb300.dat',STATUS='OLD') READ(3,10) ID,MAXS 10 FORMAT(A4,I5) IF(ID.NE.' RYB') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RYB') STOP ENDIF READ(3,30)RYBCOF 30 FORMAT(D22.16) CLOSE(3) LOADED=.TRUE. ENDIF IF(DLAMDA.GT.1D0) THEN I=3-IALPHA J=3-IBETA XL=1D0/DLAMDA ISIGN=-1 ELSE I=IALPHA J=IBETA XL=DLAMDA ISIGN=1 ENDIF XLP1=XL+1D0 C C These Gx are the functions g from Jeffrey & Onishi eqn 5.6. IF(I.EQ.1) THEN G2=-XL*(4D0+XL)/(5D0*XLP1**2) G3=-(32D0-33D0*XL+83D0*XL**2+43D0*XL**3)/(250D0*XLP1**2) ISTEP=1 ELSE G2=-(4D0*XL+1D0)/(5D0*XLP1**2) G3=-(43D0+83D0*XL-33D0*XL**2+32D0*XL**3)/(250D0*XL*XLP1**2) ISTEP=-1 ENDIF T=2D0/S EPS=1D0-4D0/(S*S) ISIGN=ISIGN*ISTEP IF(I.NE.J) THEN ZT=T*T RYB=-(G2+G3*EPS*S*S/4D0)*DLOG(EPS)-G3 MMIN=2 MMAX=MAXS ELSE ZT=T RYB=(G2+G3*EPS*S*S/4D0)*DLOG((S+2D0)/(S-2D0)) RYB=RYB-G3*S MMIN=1 MMAX=MAXS-1 ENDIF T=T*T 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 ZL=XLP1**(-M) COEF=0D0 DO 50 IQ=0,M COEF=COEF+RYBCOF(NM)*ZL NM = NM + ISTEP ZL=ZL*XL 50 CONTINUE RYB=RYB+ZT*(COEF-((G2SUM-G3SUM/DBLE(M+2))/DBLE(M))) ZT=ZT*T 100 CONTINUE RYB=RYB*DBLE(ISIGN) IF(I.EQ.J)RETURN RYB=-RYB*4D0/XLP1**2 IF(I.EQ.2)RYB=RYB*XL**2 RETURN END