Commit d3d70a24 authored by Vesa Oikonen's avatar Vesa Oikonen
Browse files

added mfEvalIntToInf()

parent 43d380af
......@@ -263,8 +263,8 @@ int mfEvalY(
/** Evaluate the function integral from 0 to specified x values.
@return Returns non-zero in case of an error.
@author Vesa Oikonen
@todo Only MF_FENGM2 and MF_FENGM2E are functional and tested.
@sa mfEvalY, mfEvalFrameY, mfEval2ndInt
@todo Only MF_FENGM2, MF_FENGM2E, erlangpdf, and exponential functions are tested.
@sa mfEvalY, mfEvalFrameY, mfEval2ndInt, mfEvalIntToInf
*/
int mfEvalInt(
/** Function code. */
......@@ -721,3 +721,79 @@ int mfEvalFrameY(
/*****************************************************************************/
/*****************************************************************************/
/** Evaluate the function integral from specified x to infinity.
@return Returns non-zero in case of an error.
@author Vesa Oikonen
@todo Only exponential functions are functional and tested.
@sa mfEvalInt, mfEvalY
*/
int mfEvalIntToInf(
/** Function code. */
const char *fid,
/** Nr of function parameters. */
const int parNr,
/** Array of function parameters. */
double *p,
/** X value; must be 0 or larger. */
double x,
/** Pointer for result integral. */
double *v,
/** Verbose level; if zero, then nothing is printed to stderr or stdout. */
const int verbose
) {
if(verbose>0) {
printf("%s('%s', %d, p[], %g, y, %d)\n", __func__, fid, parNr, x, verbose);
fflush(stdout);
}
if(parNr<1 || p==NULL || v==NULL) return(1);
if(!(x>=0.0)) return(2);
*v=0.0;
/* Exponential functions */
if(strcasecmp(fid, "1exp")==0) {
if(parNr<2) return(2);
double A=p[0], k=p[1];
*v=(-A/k)*exp(k*x);
return(0);
}
if(strcasecmp(fid, "2exp")==0) {
if(parNr<4) return(2);
for(int n=0; n<=1; n++) {
double A=p[n*2], k=p[n*2+1];
*v+=(-A/k)*exp(k*x);
}
return(0);
}
if(strcasecmp(fid, "3exp")==0) {
if(parNr<6) return(2);
for(int n=0; n<=2; n++) {
double A=p[n*2], k=p[n*2+1];
*v+=(-A/k)*exp(k*x);
}
return(0);
}
if(strcasecmp(fid, "4exp")==0) {
if(parNr<8) return(2);
for(int n=0; n<=3; n++) {
double A=p[n*2], k=p[n*2+1];
*v+=(-A/k)*exp(k*x);
}
return(0);
}
if(strcasecmp(fid, "5exp")==0) {
if(parNr<10) return(2);
for(int n=0; n<=4; n++) {
double A=p[n*2], k=p[n*2+1];
*v+=(-A/k)*exp(k*x);
}
return(0);
}
if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
return(10);
}
/*****************************************************************************/
/*****************************************************************************/
......@@ -33,9 +33,9 @@ static char *info[] = {
/*****************************************************************************/
/** Run unit tests to the library functions
* @author Vesa Oikonen
* @return 0 if all tests pass, otherwise >0.
* */
@author Vesa Oikonen
@return 0 if all tests pass, otherwise >0.
*/
int main(
/** Nr of arguments */
int argc,
......@@ -91,6 +91,11 @@ int main(
statusPrint(stderr, &status); statusFree(&status);
return(i);
}
i++; if((ret=test_mfEvalIntToInf(&status))!=0) {
fprintf(stderr, "failed (%d).\n", ret);
statusPrint(stderr, &status); statusFree(&status);
return(i);
}
/* rgamma */
i++; if((ret=test_igam(&status))!=0) {
fprintf(stderr, "failed (%d).\n", ret);
......
......@@ -620,3 +620,79 @@ int test_mfEvalFrameY(
/*****************************************************************************/
/*****************************************************************************/
int test_mfEvalIntToInf(
TPCSTATUS *status
) {
int verbose=0; if(status!=NULL) verbose=status->verbose;
statusSet(status, __func__, __FILE__, __LINE__, 0);
if(verbose>0) {
printf("\n=====================================\n");
printf("\n%s\n", __func__);
printf("\n=====================================\n");
}
int ret;
if(verbose>1) printf("\n testing with empty input \n");
ret=mfEvalIntToInf(NULL, 0, NULL, 0.0, NULL, verbose-2);
if(verbose>2) printf(" good, did not crash \n");
if(ret==0) return(1);
ret=mfEvalIntToInf("1exp", 2, NULL, 3.0, NULL, verbose-2);
if(verbose>2) printf(" good, did not crash \n");
if(ret==0) return(2);
{
char *fid="1exp";
if(verbose>1) printf("\n testing '%s' \n", fid);
double p1[]={100., -0.8};
double y;
double x=0.0;
double cy=125.;
ret=mfEvalIntToInf("1exp", 2, p1, x, &y, verbose-2);
if(ret!=0) return(101);
if(verbose>3) printf(" %g : %g vs %g \n", x, cy, y);
if(!doubleMatch(cy, y, 1.0E-05)) return(102);
if(verbose>2) printf("\n ok \n");
x=10.0;
cy=0.041932828;
ret=mfEvalIntToInf("1exp", 2, p1, x, &y, verbose-2);
if(ret!=0) return(103);
if(verbose>3) printf(" %g : %g vs %g \n", x, cy, y);
if(!doubleMatch(cy, y, 1.0E-05)) return(104);
if(verbose>2) printf("\n ok \n");
}
{
char *fid="5exp";
if(verbose>1) printf("\n testing '%s' \n", fid);
double p1[]={100., -0.8, 50., -0.2, 20., -0.05, 2., -0.01, 1., -0.001};
double y;
double x=0.0;
double cy=1975.;
ret=mfEvalIntToInf("5exp", 10, p1, x, &y, verbose-2);
if(ret!=0) return(111);
if(verbose>3) printf(" %g : %g vs %g \n", x, cy, y);
if(!doubleMatch(cy, y, 1.0E-05)) return(112);
if(verbose>2) printf("\n ok \n");
x=10.0;
cy=1447.505335;
ret=mfEvalIntToInf("5exp", 10, p1, x, &y, verbose-2);
if(ret!=0) return(113);
if(verbose>3) printf(" %g : %g vs %g \n", x, cy, y);
if(!doubleMatch(cy, y, 1.0E-05)) return(114);
if(verbose>2) printf("\n ok \n");
}
statusSet(status, __func__, __FILE__, __LINE__, 0);
return(0);
}
/*****************************************************************************/
/*****************************************************************************/
......@@ -20,6 +20,7 @@ extern int test_mfEvalY(TPCSTATUS *status);
extern int test_mfEvalInt(TPCSTATUS *status);
extern int test_mfEval2ndInt(TPCSTATUS *status);
extern int test_mfEvalFrameY(TPCSTATUS *status);
extern int test_mfEvalIntToInf(TPCSTATUS *status);
/*****************************************************************************/
/*****************************************************************************/
......
......@@ -40,6 +40,11 @@ extern int mfEvalFrameY(
const int sampleNr, double *x1, double *x2, double *y,
const int verbose
);
extern int mfEvalIntToInf(
const char *fid, const int parNr, double *p,
double x, double *v,
const int verbose
);
/*****************************************************************************/
/*****************************************************************************/
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment