/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:56 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dpsi.h" #include #include /* COMMON translations */ struct t_dpsic { double tol, xerr, psierr, round; long int ierr, msgoff; } dpsic; /* end of COMMON translations */ double /*FUNCTION*/ dpsi( double x) { long int _l0, i, n; double aux, big, dpsi_v, y; static LOGICAL32 first = TRUE; static double 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 double 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 double pi = 3.14159265358979323846264338327950e0; static long ntpsi = 0; static long ntapsi = 0; static double xbig = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Apsics = &apsics[0] - 1; double *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 DPSI Krogh Moved external statement up for mangle. *>> 1996-04-27 DPSI Krogh Changes to use .C. and C%%. *>> 1994-11-28 DPSI CLL Edit long data stmt for conversion to C. *>> 1994-11-18 DPSI WV Snyder get rid of block data *>> 1994-11-02 DPSI Krogh Changes to use M77CON *>> 1994-09-01 DPSI CLL Declare all variables. *>> 1994-05-19 DPSI WV Snyder JPL Make SP like DP using CHGTYP *>> 1993-07-23 DPSI WV Snyder JPL Adaptation from CMLIB *--D 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.00000D+00 * WITH WEIGHTED ERROR 5.79D-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.00000D-02 * WITH WEIGHTED ERROR 7.75D-33 * LOG WEIGHTED ERROR 32.11 * SIGNIFICANT FIGURES REQUIRED 28.88 * DECIMAL PLACES REQUIRED 32.71 * */ if (first) { dpsik( 0.0e0, -1.0e0, 0 ); first = FALSE; } if (ntpsi == 0) { big = DBL_EPSILON/FLT_RADIX; dinits( psics, 42, 0.1e0*big, &ntpsi ); dinits( apsics, 16, 0.1e0*big, &ntapsi ); xbig = 1.0e0/sqrt( big ); } dpsic.ierr = 0; y = fabs( x ); if (y <= 10.0e0) { /* DPSI(X) FOR ABS(X) .LE. 10 * */ if (x == 0.0e0) { dpsic.ierr = 1; ermsg( "DPSI", dpsic.ierr, 2 + dpsic.msgoff, "X is 0.0" , '.' ); dpsi_v = 0.0e0; dpsic.psierr = -1.0e0; return( dpsi_v ); } n = x; if (x < 0.0e0) { if (x == n) goto L_100; n -= 1; } y = x - n; n -= 1; dpsi_v = dcsevl( y + y - 1.0e0, psics, ntpsi ); if (n == 0) { dpsic.psierr = dpsic.round; return( dpsi_v ); } if (n <= 0) { n = -n; for (i = 0; i <= (n - 1); i++) { dpsi_v += -1.0e0/(x + i); } } else { /* DPSI(X) FOR X .GE. 2.0 AND X .LE. 10.0 * */ for (i = 1; i <= n; i++) { dpsi_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/tan( pi*x ); goto L_50; } } dpsic.psierr = dpsic.round; return( dpsi_v ); } /* DPSI(X) FOR ABS(X) .GT. 10.0 * */ if (x < 0.0e0) { n = x; if (x == n) goto L_100; } if (y < xbig) { aux = dcsevl( 2.0e0*powi(10.0e0/y,2) - 1.0e0, apsics, ntapsi ); } else { aux = 0.0e0; } if (x > 0.0e0) { dpsi_v = log( x ) - 0.5e0/x + aux; dpsic.psierr = dpsic.round; return( dpsi_v ); } n = x + x; if (x != n) { /* X is not a half-integer (X can't be integer here). */ y = pi/tan( pi*x ); } else { y = 0.0e0; } dpsi_v = log( fabs( x ) ) - 0.5e0/x + aux - y; /* Error estimate when Psi might be going to infinity or zero */ L_50: dpsic.psierr = fmax( dpsic.round, fabs( x*dpsic.xerr*y*y/dpsi_v ) ); if (dpsic.psierr > dpsic.tol) goto L_200; return( dpsi_v ); L_100: dpsic.ierr = 2; derm1( "DPSI", dpsic.ierr, 2 + dpsic.msgoff, "X is a negative integer" , "X", x, '.' ); dpsi_v = 0.0e0; dpsic.psierr = -1.0e0; return( dpsi_v ); L_200: dpsic.ierr = -3; ermsg( "DPSI", dpsic.ierr, 1 + dpsic.msgoff, "Desired precision not achieved because X is too near a negative integer" , ',' ); ermor( "or a zero of PSI(x)", ',' ); derv1( "TOL", dpsic.tol, ',' ); derv1( "ERR", dpsic.psierr, ',' ); derv1( "X", x, '.' ); return( dpsi_v ); } /* end of function */ void /*FUNCTION*/ dpsik( double toli, double xerri, long msgofi) { long int _l0; static LOGICAL32 first = TRUE; if (first) { dpsib(); dpsic.round = DBL_EPSILON; dpsic.tol = sqrt( dpsic.round ); first = FALSE; } if (toli <= 0.0e0) { dpsic.tol = sqrt( dpsic.round ); } else { dpsic.tol = fmax( toli, dpsic.round ); } if (xerri < 0.0e0) { dpsic.xerr = dpsic.round; } else { dpsic.xerr = xerri; } dpsic.msgoff = msgofi; return; } /* end of function */ void /*FUNCTION*/ dpsie( double *err, long *ierflg) { static LOGICAL32 first = TRUE; if (first) { dpsib(); first = FALSE; } *err = dpsic.psierr; *ierflg = dpsic.ierr; return; } /* end of function */ void /*FUNCTION*/ dpsib() { static LOGICAL32 first = TRUE; if (first) { first = FALSE; dpsic.ierr = 0; dpsic.psierr = 0.0; dpsic.tol = -1.0e0; } return; } /* end of function */