1 REAL FUNCTION SSINPX (X) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 1996-03-30 SSINPX Krogh Added external statement. 6C>> 1996-01-29 SCOSPX WVS JPL Add better acknowledgement for origins. 7C>> 1994-10-20 SSINPX Krogh Changes to use M77CON 8C>> 1994-04-22 SSINPX CLL Made SP and DP codes similar. 9C>> 1993-05-06 SSINPX WVS JPL Convert from NSWC to Math 77 10c--S replaces "?": ?SINPX, ?ERM1 11C ---------------------------------------------------------------------- 12c This procedure was originally procedure SIN1 from the Naval Surface 13c Warfare Center library. 14C ---------------------------------------------------------------------- 15C 16C EVALUATION OF SIN(PI*X) 17C 18C -------------- 19C 20C THE EXPANSION FOR SIN(PI*A) (ABS(A) .LE. 1/4) USING A1,...,A13 21C IS ACCURATE TO WITHIN 2 UNITS OF THE 40-TH SIGNIFICANT DIGIT, AND 22C THE EXPANSION FOR COS(PI*A) (ABS(A) .LE. 1/4) USING B1,...,B13 23C IS ACCURATE TO WITHIN 4 UNITS OF THE 40-TH SIGNIFICANT DIGIT. 24c 25c The polynomials using coefficients SA0, ... SA6 and SB1, ..., SB6 26c give approximations whose largest observed relative error in the 27c relevant intervals is 0.129e-14. 28c We will use this latter approximation when the machine epsilon 29c is larger than 0.2e-14. 30c ---------------------------------------------------------------------- 31 external R1MACH 32 integer N 33 real R1MACH 34 REAL A, BIG, CUTOFF, EPS, PI, T, W, X 35 REAL A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, 36 * A11, A12, A13 37 REAL B1, B2, B3, B4, B5, B6, B7, B8, B9, B10, 38 * B11, B12, B13 39 REAL SA0, SA1, SA2, SA3, SA4, SA5, SA6 40 REAL SB1, SB2, SB3, SB4, SB5, SB6 41 SAVE BIG, EPS 42c ----------------------- 43 PARAMETER ( PI = 3.141592653589793238462643383279502884197E+00 ) 44 parameter ( CUTOFF = 0.2e-14 ) 45c ----------------------- 46 data EPS / -1.0e0 / 47 DATA SA0 /.314159265358979E+01/, SA1 /-.516771278004995E+01/, 48 * SA2 /.255016403987327E+01/, SA3 /-.599264528932149E+00/, 49 * SA4 /.821458689493251E-01/, SA5 /-.737001831310553E-02/, 50 * SA6 /.461514425296398E-03/ 51 DATA SB1 /-.493480220054460E+01/, SB2 /.405871212639605E+01/, 52 * SB3 /-.133526276691575E+01/, SB4 /.235330543508553E+00/, 53 * SB5 /-.258048861575714E-01/, SB6 /.190653140279462E-02/ 54 DATA A1 /-.1028083791780141522795259479153765743002E+00/, 55 * A2 / .3170868848763100170457042079710451905600E-02/, 56 * A3 /-.4657026956105571623449026167864697920000E-04/, 57 * A4 / .3989844942879455643410226655783424000000E-06/, 58 * A5 /-.2237397227721999776371894030796800000000E-08/, 59 * A6 / .8847045483056962709715066675200000000000E-11/, 60 * A7 /-.2598715447506450292885585920000000000000E-13/, 61 * A8 / .5893449774331011070033920000000000000000E-16/, 62 * A9 /-.1062975472045522550784000000000000000000E-18/, 63 * A10 / .1561182648301780992000000000000000000000E-21/, 64 * A11 /-.1903193516670976000000000000000000000000E-24/, 65 * A12 / .1956617650176000000000000000000000000000E-27/, 66 * A13 /-.1711276032000000000000000000000000000000E-30/ 67c ----------------------- 68 DATA B1 /-.3084251375340424568385778437461297229882E+00/, 69 * B2 / .1585434424381550085228521039855226435920E-01/, 70 * B3 /-.3259918869273900136414318317506279360000E-03/, 71 * B4 / .3590860448591510079069203991239232000000E-05/, 72 * B5 /-.2461136950494199754009084061808640000000E-07/, 73 * B6 / .1150115912797405152263195572224000000000E-09/, 74 * B7 /-.3898073171259675439899172864000000000000E-12/, 75 * B8 / .1001886461636271969091584000000000000000E-14/, 76 * B9 /-.2019653396886572027084800000000000000000E-17/, 77 * B10 / .3278483561466560512000000000000000000000E-20/, 78 * B11 /-.4377345082051788800000000000000000000000E-23/, 79 * B12 / .4891532381388800000000000000000000000000E-26/, 80 * B13 /-.4617089843200000000000000000000000000000E-29/ 81c ----------------------- 82 IF (eps .LT. 0.0e0) then 83 eps = R1MACH(3) 84 BIG = 1.0e0/eps 85 endif 86c ----------------------- 87 A = ABS(X) 88 IF (A .GE. BIG) THEN 89 CALL SERM1 ('SSINPX',1,2, 90 1 'No precision because ABS(X) is too large','X',X,'.') 91 SSINPX = 0.E0 92 RETURN 93 END IF 94 N = A 95 T = N 96 A = A - T 97 IF (A .GT. 0.75E0) GO TO 20 98 IF (A .LT. 0.25E0) GO TO 21 99C 100C 0.25 .LE. A .LE. 0.75 101C 102 A = 0.25E0 + (0.25E0 - A) 103 if(eps .lt. cutoff) then 104 T = 16.E0*A*A 105 SSINPX = (((((((((((((B13*T + B12)*T + B11)*T + B10)*T + B9)*T+ 106 * B8)*T + B7)*T + B6)*T + B5)*T + B4)*T + B3)*T + 107 * B2)*T + B1)*T + 0.5E0) + 0.5E0 108 else 109 T = A*A 110 SSINPX = ((((((SB6*T + SB5)*T + SB4)*T + SB3)*T + SB2)*T 111 * + SB1)*T + 0.5E0) + 0.5E0 112 endif 113 GO TO 30 114C 115C A .LT. 0.25 OR A .GT. 0.75 116C 117 20 A = 0.25E0 + (0.75E0 - A) 118 21 continue 119 if(eps .lt. cutoff) then 120 T = 16.E0*A*A 121 W = (((((((((((((A13*T + A12)*T + A11)*T + A10)*T + A9)*T + 122 * A8)*T + A7)*T + A6)*T + A5)*T + A4)*T + A3)*T + 123 * A2)*T + A1)*T + 0.5E0) + 0.5E0 124 SSINPX = PI*A*W 125 else 126 T = A*A 127 SSINPX = ((((((SA6*T + SA5)*T + SA4)*T + SA3)*T + SA2)*T 128 * + SA1)*T + SA0)*A 129 endif 130C 131C TERMINATION 132C 133 30 IF (X .LT. 0.0E0) SSINPX = - SSINPX 134 IF (MOD(N,2) .NE. 0) SSINPX = - SSINPX 135 RETURN 136 END 137