INTERNAL: Elliptic integral of first kind
Anonym
Update 1/21/2013: This might be a good example of a FORTRAN program. Sometimes people come in asking about FORTRAN and an example might be helpful. Leaving internal.
Topic:
Elliptic integral of first kind
Discussion:
Fortran example of Elliptic integral of first kind
Solution:
C
C ..................................................................
C
C FUNCTION CEL1
C
C PURPOSE
C CALCULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND
C
C USAGE
C RES = CEL1(AK)
C
C DESCRIPTION OF PARAMETERS
C RES - RESULT VALUE
C AK - MODULUS (INPUT)
C
C REMARKS
C THE RESULT IS SET TO 1.E38 IF ABS(AK) GE 1
C FOR MODULUS AK AND COMPLEMENTARY MODULUS CK,
C EQUATION AK*AK+CK*CK=1.0 IS USED.
C AK MUST BE IN THE RANGE -1 TO +1
C
C
C METHOD
C DEFINITION
C CEL1(AK)=INTEGRAL(1/SQRT((1+T*T)*(1+(CK*T)**2)), SUMMED
C OVER T FROM 0 TO INFINITY).
C EQUIVALENT ARE THE DEFINITIONS
C CEL1(AK)=INTEGRAL(1/(COS(T)SQRT(1+(CK*TAN(T))**2)),SUMMED
C OVER T FROM 0 TO PI/2),
C CEL1(AK)=INTEGRAL(1/SQRT(1-(AK*SIN(T))**2),SUMMED OVER T
C FROM 0 TO PI/2), WHERE K=SQRT(1.-CK*CK).
C EVALUATION
C LANDENS TRANSFORMATION IS USED FOR CALCULATION.
C REFERENCE
C R.BULIRSCH, 'NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS
C AND ELLIPTIC FUNCTIONS', HANDBOOK SERIES SPECIAL FUNCTIONS,
C NUMERISCHE MATHEMATIK VOL. 7, 1965, PP. 78-90.
C
C ..................................................................
C
FUNCTION CEL1(AK)
C
ARI=2.
GEO=(0.5-AK)+0.5
GEO=GEO+GEO*AK
RES=0.5
IF(GEO)1,2,4
1 IER=1
2 RES=1.E38
GO TO 100
C
3 GEO=GEO*AARI
4 GEO=SQRT(GEO)
GEO=GEO+GEO
AARI=ARI
ARI=ARI+GEO
RES=RES+RES
IF(GEO/AARI-0.9999)3,5,5
5 RES=RES/ARI*6.283185E0
C
100 CEL1 = RES
RETURN
END
C
C ..................................................................
C
C FUNCTION EXPI
C
C PURPOSE
C COMPUTES THE EXPONENTIAL INTEGRAL -EI(-X)
C
C USAGE
C RES = EXPI(X)
C
C DESCRIPTION OF PARAMETERS
C X - ARGUMENT OF EXPONENTIAL INTEGRAL
C RES - RESULT VALUE
C
C REMARKS
C X GT 87 MAY CAUSE UNDERFLOW
C WITH THE EXPONENTIAL FUNCTION
C FOR X = 0 THE RESULT VALUE IS SET TO -1.E37
C
C METHOD
C DEFINITION
C RES=INTEGRAL(EXP(-T)/T, SUMMED OVER T FROM X TO INFINITY).
C EVALUATION
C THREE DIFFERENT RATIONAL APPROXIMATIONS ARE USED IN THE
C RANGES 1 LE X, X LE -9 AND -9 LT X LE -3 RESPECTIVELY,
C A POLYNOMIAL APPROXIMATION IS USED IN -3 LT X LT 1.
C
C ..................................................................
C
FUNCTION EXPI(X)
C
IF(X-1.)2,1,1
1 Y=1./X
AUX=1.-Y*(((Y+3.377358E0)*Y+2.052156E0)*Y+2.709479E-1)/((((Y*
11.072553E0+5.716943E0)*Y+6.945239E0)*Y+2.593888E0)*Y+2.709496E-1)
RES=AUX*Y*EXP(-X)
GO TO 100
C
2 IF (X+3.) 6,6,3
C
3 AUX=(((((((7.122452E-7*X-1.766345E-6)*X+2.928433E-5)*X-2.335379E-4
1)*X+1.664156E-3)*X-1.041576E-2)*X+5.555682E-2)*X-2.500001E-1)*X
2+9.999999E-1
RES=-1.E37
IF(X) 4,100,4
4 RES=X*AUX-ALOG(ABS(X))-5.772157E-1
GO TO 100
C
6 IF (X+9.)8,8,7
C
7 AUX=1.-((((5.176245E-2*X+3.061037E0)*X+3.243665E1)*X+2.244234E2)*X
1+2.486697E2)/((((X+3.995161E0)*X+3.893944E1)*X+2.263818E1)*X
2+1.807837E2)
GOTO 9
8 Y=9./X
AUX=1.-Y*(((Y+7.659824E-1)*Y-7.271015E-1)*Y-1.080693E0)/((((Y
1*2.518750E0+1.122927E1)*Y+5.921405E0)*Y-8.666702E0)*Y-9.724216E0)
9 RES=AUX*EXP(-X)/X
C
100 EXPI = RES
RETURN
END
C
C ..................................................................
C
C FUNCTION SINI
C
C PURPOSE
C COMPUTES THE SINE INTEGRAL
C
C USAGE
C SI = SINI(X)
C
C DESCRIPTION OF PARAMETERS
C SI - THE RESULTANT VALUE SI(X)
C X - THE ARGUMENT OF SI(X) AND CI(X)
C
C REMARKS
C THE ARGUMENT VALUE REMAINS UNCHANGED
C
C METHOD
C DEFINITION
C SI(X)=INTEGRAL(SIN(T)/T)
C EVALUATION
C REDUCTION OF RANGE USING SYMMETRY.
C DIFFERENT APPROXIMATIONS ARE USED FOR ABS(X) GREATER
C THAN 4 AND FOR ABS(X) LESS THAN 4.
C REFERENCE
C LUKE AND WIMP, 'POLYNOMIAL APPROXIMATIONS TO INTEGRAL
C TRANSFORMS', MATHEMATICAL TABLES AND OTHER AIDS TO
C COMPUTATION, VOL. 15, 1961, ISSUE 74, PP. 174-178.
C
C ..................................................................
C
C ..................................................................
C
C FUNCTION COSI
C
C PURPOSE
C COMPUTES THE COSINE INTEGRAL
C
C USAGE
C CI = COSI(X)
C
C DESCRIPTION OF PARAMETERS
C CI - THE RESULTANT VALUE CI(X)
C X - THE ARGUMENT OF SI(X) AND CI(X)
C
C REMARKS
C THE ARGUMENT VALUE REMAINS UNCHANGED
C
C METHOD
C DEFINITION
C CI(X)=INTEGRAL(COS(T)/T)
C EVALUATION
C REDUCTION OF RANGE USING SYMMETRY.
C DIFFERENT APPROXIMATIONS ARE USED FOR ABS(X) GREATER
C THAN 4 AND FOR ABS(X) LESS THAN 4.
C REFERENCE
C LUKE AND WIMP, 'POLYNOMIAL APPROXIMATIONS TO INTEGRAL
C TRANSFORMS', MATHEMATICAL TABLES AND OTHER AIDS TO
C COMPUTATION, VOL. 15, 1961, ISSUE 74, PP. 174-178.
C
C ..................................................................
C
C
FUNCTION SINI(X)
parameter PI = 3.141592653
isini = .true.
goto 1000
c
entry cosi(x)
isini = .false.
C
1000 Z=ABS(X)
IF(Z-4.)1,1,4
1 Y=(4.-Z)*(4.+Z)
SI=-1.570797E0
C
IF(Z)3,2,3
2 CI=-1.E37
GO TO 100
C
3 SI=X*(((((1.753141E-9*Y+1.568988E-7)*Y+1.374168E-5)*Y+6.939889E-4)
1*Y+1.964882E-2)*Y+4.395509E-1+SI/X)
CI=((5.772156E-1+ALOG(Z))/Z-Z*(((((1.386985E-10*Y+1.584996E-8)*Y
1+1.725752E-6)*Y+1.185999E-4)*Y+4.990920E-3)*Y+1.315308E-1))*Z
GO TO 100
C
4 SI=SIN(Z)
Y=COS(Z)
Z=4./Z
U=((((((((4.048069E-3*Z-2.279143E-2)*Z+5.515070E-2)*Z-7.261642E-2)
1*Z+4.987716E-2)*Z-3.332519E-3)*Z-2.314617E-2)*Z-1.134958E-5)*Z
2+6.250011E-2)*Z+2.583989E-10
V=(((((-5.108699E-3*Z+2.819179E-2)*Z-6.537283E-2)*Z
1+7.902034E-2)*Z-4.400416E-2)*Z-7.945556E-3)*Z+2.601293E-2
V=(((V*Z-3.764000E-4)*Z-3.122418E-2)*Z-6.646441E-7)*Z+2.500000E-1
CI=Z*(SI*V-Y*U)
SI=-Z*(SI*U+Y*V)
IF (X) 5,100,100
5 SI=pi-SI
C
100 if (isini) then
sini = si+pi/2.
else
sini = ci
end if
RETURN
END
FUNCTION VOIGT(Y,X)
C
C
C...SEE J.Q.S.R.T. VOLUME 7, PAGE 85.
C
C
DIMENSION W(10),T(10),C(34)
C
DATA C/
1 0.199999999997222E0,-0.184000000003000E0, 0.155839999996502E0,
2 -0.121664000004399E0, 0.087708159994039E0,-0.058514124808691E0,
3 0.036215730162391E0,-0.020849765439804E0, 0.011196011634627E0,
4 -0.56231896167109E-2, 0.26487634172265E-2,-0.11732670757704E-2,
5 0.4899519978088E-3, -0.1933630801528E-3, 0.722877446788E-4,
6 -0.256555124979E-4, 0.86620736841E-5, -0.27876379719E-5,
7 0.8566873627E-6, -0.2518433784E-6, 0.709360221E-7,
8 -0.191732257E-7, 0.49801256E-8, -0.12447734E-8,
9 0.2997777E-9, -0.696450E-10, 0.156262E-10,
1 -0.33897E-11, 0.7116E-12, -0.1447E-12,
1 0.285E-13, -0.55E-14, 0.10E-14,
2 -0.2E-15/
C
DATA W/
1 4.62243670E-1, 2.86675505E-1,1.09017206E-1,2.48105209E-2,
2 3.24377334E-3, 2.28338636E-4,7.80255648E-6,1.08606937E-7,
3 4.39934099E-10,2.22939365E-13/
C
DATA T/
1 0.245340708,0.737473729,1.23407622,1.73853771,2.25497400,
2 2.78880606, 3.34785457, 3.94476404,4.60368245,5.38748089/
C
C
F3(Z) = EXP(Z*Z-X*X)
C
Y2 = Y*Y
IF(Y.LT.1.0.AND.X.LT.4.0.OR.Y.LT.1.8/(X+1.0)) GO TO 300
IF(Y.LT.2.5.AND.X.LT.4.0) GO TO 200
C
100 G = 0.0E0
DO 101 I=1,10
G=G+(1.0E0/((X-T(I))**2+Y2)+1.0E0/((X+T(I))**2+Y2))*W(I)
101 CONTINUE
VOIGT = 0.318309886*Y*G
RETURN
C
C
200 G = 0.0E0
DO 201 I=1,10
R = T(I)-X
S = T(I)+X
G=G+(4.*T(I)**2-2.)*(R* ATAN(R/Y)+S* ATAN(S/Y)-.5*Y*(ALOG(Y2+R*R)+
1 ALOG(Y2+S*S)))*W(I)
201 CONTINUE
VOIGT = 0.318309886*G
RETURN
C
C
300 IF((X*X-Y2).GT.100.E0) GO TO 2
U1 = EXP(-X*X+Y2)*COS(2.*X*Y)
GO TO 5
2 U1 = 0.0E0
5 IF(X.GT.5.0) GO TO 1000
C
C
BN01 = 0.0E0
BN02 = 0.0E0
X1 = X/5.0E0
COEF = 4.0E0*X1*X1-2.0E0
DO 20 I=1,34
II = 35-I
BN = COEF*BN01-BN02+C(II)
BN02 = BN01
BN01 = BN
20 CONTINUE
30 F = X1*(BN-BN02)
40 DN01 = 1.E0 - 2.E0*X*F
1100 DN02 = F
GO TO 1200
C
1000 Z = 1.E0/(2.E0*X*X)
DN01=-Z*(1.E0+3.E0*Z*(1.E0+5.E0*Z*(1.E0+7.E0*Z*(1.E0+9.E0*Z*
1 (1.E0+11.E0*Z*(1.E0+13.E0*Z))))))
DN02 = (1.E0-DN01)/(2.E0*X)
1200 FUNCT = Y*DN01
IF(Y.LE.1.0E-8) GO TO 2500
Q = 1.0
YN = Y
DO 2000 I=2,50
DN = (X*DN01+DN02)*(-2.)/ FLOAT(I)
DN02 = DN01
DN01 = DN
IF(MOD(I,2)) 2000,2000,1500
1500 Q = -Q
YN = YN*Y2
G = DN*YN
FUNCT = FUNCT + Q*G
IF( ABS(G/FUNCT).LE.1.0E-8) GO TO 2500
2000 CONTINUE
2500 VOIGT = U1 - 1.12837917*FUNCT
RETURN
END