DECLARE SUB CISIA (X#, CI#, SI#) '****************************************************************** '* Purpose: This program computes the cosine and sine * '* integrals using subroutine CISIA * '* Input : x --- Argument of Ci(x) and Si(x) * '* Output: CI --- Ci(x) * '* SI --- Si(x) * '* Example: * '* x Ci(x) Si(x) * '* ------------------------------------ * '* 0.0 - i .00000000 * '* 5.0 -.19002975 1.54993124 * '* 10.0 -.04545643 1.65834759 * '* 20.0 .04441982 1.54824170 * '* 30.0 -.03303242 1.56675654 * '* 40.0 .01902001 1.58698512 * '* -------------------------------------------------------------- * '* REFERENCE: "Fortran Routines for Computation of Special * '* Functions jin.ece.uiuc.edu/routines/routines.html" * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '****************************************************************** DEFDBL A-H, O-Z DEFINT I-N ' PROGRAM MCISIA PRINT INPUT " Please enter x: ", X PRINT " x Ci(x) Si(x) " PRINT "------------------------------------" CALL CISIA(X, CI, SI) IF X <> 0! THEN PRINT " "; X; " "; CI; " "; SI IF X = 0! THEN PRINT " "; X; " i"; " O" END 'end of file mcisia.bas SUB CISIA (X, CI, SI) ' ============================================= ' Purpose: Compute cosine and sine integrals ' Si(x) and Ci(x) ( x ò 0 ) ' Input : x --- Argument of Ci(x) and Si(x) ' Output: CI --- Ci(x) ' SI --- Si(x) ' ============================================= DIM BJ(100) P2 = 1.570796326794897# EL = .5772156649015329# EPS = .000000000000001# X2 = X * X IF (X = 0!) THEN CI = -1000000000000000# SI = 0! ELSEIF (X <= 16!) THEN XR = -.25 * X2 CI = EL + LOG(X) + XR FOR K = 2 TO 40 XR = -.5 * XR * (K - 1) / (K * K * (2 * K - 1)) * X2 CI = CI + XR IF (ABS(XR) < ABS(CI) * EPS) THEN GOTO 15 NEXT K 15 XR = X SI = X FOR K = 1 TO 40 XR = -.5# * XR * (2 * K - 1) / K / (4 * K * K + 4 * K + 1) * X2 SI = SI + XR IF (ABS(XR) < ABS(SI) * EPS) THEN GOTO 20'return NEXT K ELSEIF (X <= 32!) THEN M = INT(47.2 + .82 * X) XA1 = 0! XA0 = .000000000000001# FOR K = M TO 1 STEP -1 XA = 4! * K * XA0 / X - XA1 BJ(K) = XA XA1 = XA0 XA0 = XA NEXT K XS = BJ(1) FOR K = 3 TO M STEP 2 XS = XS + 2! * BJ(K) NEXT K BJ(1) = BJ(1) / XS FOR K = 2 TO M BJ(K) = BJ(K) / XS NEXT K XR = 1! XG1 = BJ(1) FOR K = 2 TO M XR = .25 * XR * (2! * K - 3!) ^ 2 / ((K - 1!) * (2! * K - 1!) ^ 2) * X XG1 = XG1 + BJ(K) * XR NEXT K XR = 1! XG2 = BJ(1) FOR K = 2 TO M XR = .25 * XR * (2! * K - 5!) ^ 2 / ((K - 1!) * (2! * K - 3!) ^ 2) * X XG2 = XG2 + BJ(K) * XR NEXT K XCS = COS(X / 2!) XSS = SIN(X / 2!) CI = EL + LOG(X) - X * XSS * XG1 + 2 * XCS * XG2 - 2 * XCS * XCS SI = X * XCS * XG1 + 2 * XSS * XG2 - SIN(X) ELSE XR = 1! XF = 1! FOR K = 1 TO 9 XR = -2! * XR * K * (2 * K - 1) / X2 XF = XF + XR NEXT K XR = 1! / X XG = XR FOR K = 1 TO 8 XR = -2! * XR * (2 * K + 1) * K / X2 XG = XG + XR NEXT K CI = XF * SIN(X) / X - XG * COS(X) / X SI = P2 - XF * COS(X) / X - XG * SIN(X) / X END IF 20 END SUB