1 REAL FUNCTION SCI (X) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5C>> 1998-10-29 SCI Krogh Moved external statement up for mangle. 6c>> 1996-03-30 SCI Krogh Added external statements. 7C>> 1995-11-22 SCI Krogh Removed multiple entry for C conversion. 8C>> 1995-11-03 SCI Krogh Removed blanks in numbers for C conversion. 9c>> 1994-10-20 SCI Krogh Changes to use M77CON 10c>> 1990-01-23 SCI CLL Using name SCIN for result in entry SCIN. 11c>> 1989-03-14 Original W. V. Snyder at JPL 12C 13C COMPUTE THE COSINE INTEGRAL OF X = 14C INTEGRAL FROM X TO INFINITY OF -(COS(T)/T) DT = 15C GAMMA + LOG(ABS(X)) + INTEGRAL FROM 0 TO X OF ((COS(T)-1)/T DT), 16C WHERE GAMMA IS EULER'S CONSTANT. 17C 18C FOR ABS(X)<16, USE A CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE 19C Z=X/16 TO EVALUATE (CI(X)-LOG(ABS(X))-GAMMA)/(Z*Z), THEN MULTIPLY 20C THE RESULT BY Z*Z AND ADD LOG(ABS(X)) AND GAMMA. LOSS OF ABOUT 21C TWO DIGITS OCCURS NEAR ABS(X)=16 DUE TO CANCELLATION. 22C 23C FOR ABS(X).GE.16, USE CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE 24C WHERE Z=16/X ARE USED TO COMPUTE F(X)/Z AND G(X)/(Z*Z). THEN 25C CI(X)=F(X)*SIN(X)-G(X)*COS(X). 26C 27C THIS ALGORITHM YIELDS AT MOST 15 DIGITS OF PRECISION. 28C 29C--S replaces "?": ?CI, ?CII, ?CIN, ?CPVAL, ?ERM1 30C 31 EXTERNAL SCII 32 REAL X, SCII 33 SCI = SCII(.true., X) 34 end 35C 36 REAL FUNCTION SCIN (X) 37 EXTERNAL SCII 38 REAL X, SCII 39 SCIN = SCII(.false., X) 40 end 41c 42 REAL FUNCTION SCII (LDCI, X) 43 LOGICAL LDCI 44 INTEGER ERRLEV 45 PARAMETER (ERRLEV=0) 46 EXTERNAL SCPVAL 47 REAL C(23), SCPVAL, F(13), G(13), GAMMA,X, Z,ZW 48 parameter(GAMMA = 0.577215664901533E0) 49 DATA C/ 50 1 + 0.5E0, + 0.5E0, 51 2 +14.992589367813409E0, -19.386124096607770E0, 52 3 +12.741870869758071E0, - 8.107903970562531E0, 53 4 + 4.862022348500627E0, - 2.497505088539025E0, 54 5 + 1.008660787358110E0, - 0.312080924825428E0, 55 6 + 0.074678255294576E0, - 0.014110865253535E0, 56 7 + 0.002152046752074E0, - 0.000270212331184E0, 57 8 + 0.000028416945498E0, - 0.000002540125611E0, 58 9 + 0.000000195437144E0, - 0.000000013084020E0, 59 A + 0.000000000769379E0, - 0.000000000040066E0, 60 B + 0.000000000001861E0, - 0.000000000000078E0, 61 C + 0.000000000000003E0/ 62 DATA F/ 63 1 + 0.5E0, + 0.5E0, 64 2 + 0.062263729028927E0, - 0.000233756041393E0, 65 3 + 0.000002453755677E0, - 0.000000058670317E0, 66 4 + 0.000000002356196E0, - 0.000000000136096E0, 67 5 + 0.000000000010308E0, - 0.000000000000964E0, 68 6 + 0.000000000000107E0, - 0.000000000000014E0, 69 7 + 0.000000000000002E0/ 70 DATA G/ 71 1 + 0.5E0, + 0.5E0, 72 2 + 0.003862856096703E0, - 0.000042644182622E0, 73 3 + 0.000000724995950E0, - 0.000000023468225E0, 74 4 + 0.000000001169202E0, - 0.000000000079604E0, 75 5 + 0.000000000006875E0, - 0.000000000000717E0, 76 6 + 0.000000000000087E0, - 0.000000000000012E0, 77 7 + 0.000000000000002E0/ 78C 79 IF (LDCI) then 80 IF (X.LE.0.0) THEN 81 CALL SERM1('SCI',1,ERRLEV,'Argument not positive','X',X,'.') 82C Provide a value in case the error message processor returns. 83 SCII = 0.0 84 ELSE IF (X.LT.16.0) THEN 85 Z = X/16.0 86 ZW = Z*Z 87 SCII = LOG(X) - ZW*SCPVAL(C,20,ZW) + GAMMA 88 ELSE 89 Z = 16.0/X 90 ZW = Z*Z 91 SCII = Z*(SCPVAL(F,10,ZW)*SIN(X) - Z*SCPVAL(G,10,ZW)*COS(X)) 92 END IF 93 ELSE 94C 95C Evaluate the entire function 96C Cin(x) = Integral from 0 to x of ((cos(t) - 1) / t) dt. 97C 98 IF (ABS(X).LT.16.0) THEN 99 Z = X/16.0 100 ZW = Z*Z 101 SCII = ZW*SCPVAL(C,20,ZW) 102 ELSE 103 Z = 16.0/X 104 ZW = Z*Z 105 SCII = LOG(ABS(X)) + GAMMA - 106 1 Z*(SCPVAL(F,10,ZW)*SIN(X) - Z*SCPVAL(G,10,ZW)*COS(X)) 107 END IF 108 END IF 109 RETURN 110 END 111