/* CDF for the S distribution. Last modified on 01/07/02 by Sergej V. Aksenov. Return codes: 0 - ok, 1-8 - CVODE error, 9 - CVODE gives negative solution, 10-15 - lerchphi() errors. */ #include static double *pars1, *f01, xgoal, scale, prec1; static double fun(double *f, int *flag) { double xcurr, fscaled; extern double *pars1, xgoal, scale; fscaled = (*f)*scale; *flag = sdqnt(pars1,f01,&fscaled,&xcurr,&prec1); if (*flag) return 1.0; else return (xgoal - xcurr); } int sdcdf(double *sdpar, double *f0, double *data, long lendat, double *cdf, double *lprec, double *ztol, double *rtol, double *atol, int *meth, int *steps) { double fasym = 0.0, ax, *leftdata, *leftres, *rightres, g, x0, *xend, fa, fb, scale1, funu, funl, zerosol, fascaled, fbscaled; int i, iax, im, status1, status2, indsol, oindsol; extern double *pars1, xgoal, scale, prec1; double zeroin_(); prec1 = *lprec; g = sdpar[0]; x0 = sdpar[3]; /* Asymptotic x; needed not to request solution of the S-d ODE past it. */ status1 = 0; if (g < 1.0) { status1 = sdqnt(sdpar,f0,&fasym,&ax,lprec); if (status1) return (status1+9); } else ax = data[0]; /* Calculate positions in data array corresponding to points greater or equal to x0 and asymptotic x: x0 index is needed to separate data array into two parts, to the right and left of x0; asymptotic x index is needed not to pass data points past (less) than ax. */ /* iax and im are preinitialized to maximum indices plus 1. */ iax = lendat; im = lendat; /* Loop from first to last index. */ for (i = 0; i <= lendat-1; i++) { /* Check for data[i] only of iax and im have not been changed since initialization. */ if (iax == lendat) { if (data[i] >= ax) iax = i; } if (im == lendat) { if (data[i] >= x0) im = i; } } /* Allocate space for portions of data and cdf between asymptotic x and x0. (index runs from iax to im-1). */ if (im-iax > 0) { leftdata = (double *) malloc((im-iax) * sizeof(double)); leftres = (double *) malloc((im-iax) * sizeof(double)); } /* Allocate space for cdf between x0 and maximum. (index runs from im to lendat-1). */ if (lendat-im > 0) rightres = (double *) malloc((lendat-im) * sizeof(double)); /* Arrange data to the left of median in order of increasing absolute value (lower index - closer to median; required by the ODE solver). */ for (i = 0; i <= im - iax - 1; i++) { leftdata[i] = data[im-i-1]; } /* Call ODE solver for portion of data to the right of x0. */ if (lendat-im > 0) { indsol = lendat-im-1; xend = &data[im]; status2 = 0; status2 = sdode(sdpar,f0,xend,&indsol,rightres,rtol,atol,meth,steps); if (status2) { if (im-iax > 0) { free(leftdata); free(leftres); } free(rightres); return status2; } } /* Call ODE solver (and solve nonlinear equation if needed) for portion of data to the left of x0. */ if (im-iax > 0) { oindsol = indsol = im-iax-1; xend = leftdata; status2 = 0; status2 = sdode(sdpar,f0,xend,&indsol,leftres,rtol,atol,meth,steps); if (status2) { /* If g>=1.0 then CVODE error has nothing to do with singularity, so exit with error. */ if (g >= 1.0) { free(leftdata); free(leftres); if (lendat-im > 0) free(rightres); return status2; } /* If g<1.0 then CVODE error might be because of singularity. So bracket the zero, then find zero of a nonlinear eq. */ else { /* In case of CVODE error, indsol contains the index of last correctly solved point. */ i = indsol; while (i < oindsol) { /* Target x where solution is desired. */ i++; xgoal = *(xend+i); /* if xgoal is the asymptotic x, then solution f is simply 0.0 */ if (xgoal == ax) *(leftres+i) = 0.0; else { /* Largest interval where solution might be (in reality, CVODE never has problems here but with much lower fa and fb). */ fa = 0.01; fb = 0.1; scale1 = 100.; /* Scale f by 1 while bracketing. */ scale = 1.; pars1 = sdpar; f01 = f0; status1 = 0; funl = fun(&fa,&status1); if (status1) { free(leftdata); free(leftres); free(rightres); return (status1+9); } funu = fun(&fb,&status1); if (status1) { free(leftdata); free(leftres); free(rightres); return (status1+9); } /* Positive funu means true solution is larger than f=0.1 and CVODE failure is not from singularity (f<<1). */ if (funu >= 0.0) { free(leftdata); free(leftres); free(rightres); return status2; } else { /* Lower f's while fun is negative, until funl becomes positive and zero is bracketed. scale1 is increased so that in the end f's are between 1 and 10 */ while (funu < 0.0 && funl < 0.0) { scale1 *= 10; fa /= 10; fb /= 10; funu = funl; status1 = 0; funl = fun(&fa,&status1); if (status1) { free(leftdata); free(leftres); free(rightres); return (status1+9); } } /* Update scale so that fun() works with scaled f's. */ scale = 1./scale1; /* Solve nonlinear eq. with scaled f's. */ fascaled = fa*scale1; fbscaled = fb*scale1; status1 = 0; zerosol = scale*zeroin_(&fascaled, &fbscaled, &fun, ztol, status1); if (status1) { free(leftdata); free(leftres); free(rightres); return (status1+9); } *(leftres+i) = zerosol; } /* else fun>=0 */ } /* else xgoal==ax */ } /* while i=1 */ } /* if status2 */ } /* if im>iax */ /* Fill in cdf array in proper order. */ for (i = 0; i <= lendat - 1; i++) { if ((0 <= i) && (i <= iax-1)) cdf[i] = 0.0; if ((i >= iax) && (i <= im-1)) cdf[i] = leftres[im-1-i]; if (i >= im) cdf[i] = rightres[i-im]; } /* Clean up. */ if (im-iax > 0) { free(leftdata); free(leftres); } if (lendat-im > 0) free(rightres); return 0; }