/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:59 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "spsi.h" #include #include /* COMMON translations */ struct t_spsic { float tol, xerr, psierr, round; long int ierr, msgoff; } spsic; /* end of COMMON translations */ float /*FUNCTION*/ spsi( float x) { long int _l0, i, n; float aux, big, spsi_v, y; static LOGICAL32 first = TRUE; static float psics[42]={-.38057080835217921520437677667039e-1, .49141539302938712748204699654277e0,-.56815747821244730242892064734081e-1, .83578212259143131362775650747862e-2,-.13332328579943425998079274172393e-2, .22031328706930824892872397979521e-3,-.37040238178456883592889086949229e-4, .62837936548549898933651418717690e-5,-.10712639085061849855283541747074e-5, .18312839465484165805731589810378e-6,-.31353509361808509869005779796885e-7, .53728087762007766260471919143615e-8,-.92116814159784275717880632624730e-9, .15798126521481822782252884032823e-9,-.27098646132380443065440589409707e-10, .46487228599096834872947319529549e-11,-.79752725638303689726504797772737e-12, .13682723857476992249251053892838e-12,-.23475156060658972717320677980719e-13, .40276307155603541107907925006281e-14,-.69102518531179037846547422974771e-15, .11856047138863349552929139525768e-15,-.20341689616261559308154210484223e-16, .34900749686463043850374232932351e-17,-.59880146934976711003011081393493e-18, .10273801628080588258398005712213e-18,-.17627049424561071368359260105386e-19, .30243228018156920457454035490133e-20,-.51889168302092313774286088874666e-21, .89027730345845713905005887487999e-22,-.15274742899426728392894971904000e-22, .26207314798962083136358318079999e-23,-.44964642738220696772598388053333e-24, .77147129596345107028919364266666e-25,-.13236354761887702968102638933333e-25, .22709994362408300091277311999999e-26,-.38964190215374115954491391999999e-27, .66851981388855302310679893333333e-28,-.11469986654920864872529919999999e-28, .19679385886541405920515413333333e-29,-.33764488189750979801907200000000e-30, .57930703193214159246677333333333e-31}; static float apsics[16]={-.832710791069290760174456932269e-3,-.416251842192739352821627121990e-3, .103431560978741291174463193961e-6,-.121468184135904152987299556365e-9, .311369431998356155521240278178e-12,-.136461337193177041776516100945e-14, .902051751315416565130837974000e-17,-.831542997421591464829933635466e-19, .101224257073907254188479482666e-20,-.156270249435622507620478933333e-22, .296542716808903896133226666666e-24,-.674686886765702163741866666666e-26, .180345311697189904213333333333e-27,-.556901618245983607466666666666e-29, .195867922607736251733333333333e-30,-.775195892523335680000000000000e-32}; static float pi = 3.14159265358979323846264338327950e0; static long ntpsi = 0; static long ntapsi = 0; static float xbig = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Apsics = &apsics[0] - 1; float *const Psics = &psics[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1998-10-29 SPSI Krogh Moved external statement up for mangle. *>> 1996-04-27 SPSI Krogh Changes to use .C. and C%%. *>> 1994-11-28 SPSI CLL Edit long data stmt for conversion to C. *>> 1994-11-18 SPSI WV Snyder get rid of block data *>> 1994-11-02 SPSI Krogh Changes to use M77CON *>> 1994-09-01 SPSI CLL Declare all variables. *>> 1994-05-19 SPSI WV Snyder JPL Make SP like DP using CHGTYP *>> 1993-07-23 SPSI WV Snyder JPL Adaptation from CMLIB *--S replaces "?": ?PSI, ?PSIC, ?CSEVL, ?ERM1, ?INITS, ?PSIK, ?ERV1 *--& ?PSIB, ?PSIE * JUNE 1977 EDITION. W. FULLERTON, C3, LOS ALAMOS SCIENTIFIC LAB. */ /* The array APSCS(1:42) contains coefficients for the * series for PSI ON THE INTERVAL 0. TO 1.00000E+00 * WITH WEIGHTED ERROR 5.79E-32 * LOG WEIGHTED ERROR 31.24 * SIGNIFICANT FIGURES REQUIRED 30.93 * DECIMAL PLACES REQUIRED 32.05 *++ Save data by elements if ~.C. */ /* The array APSICS(1:16) contains coefficients for the * series for APSI ON THE INTERVAL 0. TO 1.00000E-02 * WITH WEIGHTED ERROR 7.75E-33 * LOG WEIGHTED ERROR 32.11 * SIGNIFICANT FIGURES REQUIRED 28.88 * DECIMAL PLACES REQUIRED 32.71 * */ if (first) { spsik( 0.0e0, -1.0e0, 0 ); first = FALSE; } if (ntpsi == 0) { big = FLT_EPSILON/FLT_RADIX; sinits( psics, 42, 0.1e0*big, &ntpsi ); sinits( apsics, 16, 0.1e0*big, &ntapsi ); xbig = 1.0e0/sqrtf( big ); } spsic.ierr = 0; y = fabsf( x ); if (y <= 10.0e0) { /* SPSI(X) FOR ABS(X) .LE. 10 * */ if (x == 0.0e0) { spsic.ierr = 1; ermsg( "SPSI", spsic.ierr, 2 + spsic.msgoff, "X is 0.0" , '.' ); spsi_v = 0.0e0; spsic.psierr = -1.0e0; return( spsi_v ); } n = x; if (x < 0.0e0) { if (x == n) goto L_100; n -= 1; } y = x - n; n -= 1; spsi_v = scsevl( y + y - 1.0e0, psics, ntpsi ); if (n == 0) { spsic.psierr = spsic.round; return( spsi_v ); } if (n <= 0) { n = -n; for (i = 0; i <= (n - 1); i++) { spsi_v += -1.0e0/(x + i); } } else { /* SPSI(X) FOR X .GE. 2.0 AND X .LE. 10.0 * */ for (i = 1; i <= n; i++) { spsi_v += 1.0e0/(y + i); } } if (x < -0.5e0) { n = x + x; if (x != n) { /* X is not a half-integer (X can't be integer here). */ y = pi/tanf( pi*x ); goto L_50; } } spsic.psierr = spsic.round; return( spsi_v ); } /* SPSI(X) FOR ABS(X) .GT. 10.0 * */ if (x < 0.0e0) { n = x; if (x == n) goto L_100; } if (y < xbig) { aux = scsevl( 2.0e0*powif(10.0e0/y,2) - 1.0e0, apsics, ntapsi ); } else { aux = 0.0e0; } if (x > 0.0e0) { spsi_v = logf( x ) - 0.5e0/x + aux; spsic.psierr = spsic.round; return( spsi_v ); } n = x + x; if (x != n) { /* X is not a half-integer (X can't be integer here). */ y = pi/tanf( pi*x ); } else { y = 0.0e0; } spsi_v = logf( fabsf( x ) ) - 0.5e0/x + aux - y; /* Error estimate when Psi might be going to infinity or zero */ L_50: spsic.psierr = fmaxf( spsic.round, fabsf( x*spsic.xerr*y*y/spsi_v ) ); if (spsic.psierr > spsic.tol) goto L_200; return( spsi_v ); L_100: spsic.ierr = 2; serm1( "SPSI", spsic.ierr, 2 + spsic.msgoff, "X is a negative integer" , "X", x, '.' ); spsi_v = 0.0e0; spsic.psierr = -1.0e0; return( spsi_v ); L_200: spsic.ierr = -3; ermsg( "SPSI", spsic.ierr, 1 + spsic.msgoff, "Desired precision not achieved because X is too near a negative integer" , ',' ); ermor( "or a zero of PSI(x)", ',' ); serv1( "TOL", spsic.tol, ',' ); serv1( "ERR", spsic.psierr, ',' ); serv1( "X", x, '.' ); return( spsi_v ); } /* end of function */ void /*FUNCTION*/ spsik( float toli, float xerri, long msgofi) { long int _l0; static LOGICAL32 first = TRUE; if (first) { spsib(); spsic.round = FLT_EPSILON; spsic.tol = sqrtf( spsic.round ); first = FALSE; } if (toli <= 0.0e0) { spsic.tol = sqrtf( spsic.round ); } else { spsic.tol = fmaxf( toli, spsic.round ); } if (xerri < 0.0e0) { spsic.xerr = spsic.round; } else { spsic.xerr = xerri; } spsic.msgoff = msgofi; return; } /* end of function */ void /*FUNCTION*/ spsie( float *err, long *ierflg) { static LOGICAL32 first = TRUE; if (first) { spsib(); first = FALSE; } *err = spsic.psierr; *ierflg = spsic.ierr; return; } /* end of function */ void /*FUNCTION*/ spsib() { static LOGICAL32 first = TRUE; if (first) { first = FALSE; spsic.ierr = 0; spsic.psierr = 0.0; spsic.tol = -1.0e0; } return; } /* end of function */