TPCCLIB
Loading...
Searching...
No Matches
Functions
tpcfunc.h File Reference

Header file for libtpcfunc. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "tpcextensions.h"
#include "tpcmodels.h"

Go to the source code of this file.

Functions

int mfEvalY (const char *fid, const int parNr, const double *p, const int sampleNr, const double *x, double *y, const int verbose)
int mfEvalInt (const char *fid, const int parNr, const double *p, const int sampleNr, const double *x, double *i, const int verbose)
int mfEval2ndInt (const char *fid, const int parNr, const double *p, const int sampleNr, const double *x, double *i, const int verbose)
int mfEvalFrameY (const char *fid, const int parNr, const double *p, const int sampleNr, const double *x1, const double *x2, double *y, const int verbose)
int mfEvalIntToInf (const char *fid, const int parNr, const double *p, const double x, double *v, const int verbose)
double igam (double a, double x)
double igamc (double a, double x)

Detailed Description

Header file for libtpcfunc.

Author
Vesa Oikonen

Definition in file tpcfunc.h.

Function Documentation

◆ igam()

double igam ( double a,
double x )
extern

Cumulative gamma distribution, or Regularized gamma function, more specifically, lower incomplete gamma function divided by gamma function.

Standard gamma distribution is assumed (Beta=1). f(a,x) = (1/Gamma(a)) * Integral(0,x)(e^-t * t^(a-1))dt

Returns
Returns the value of regularized gamma function, or NaN in case of an error.
See also
igamc, inverfc, factorial
Parameters
aShape parameter alpha; must be > 0.
xIntegral from 0 to x; must be >= 0.

Definition at line 29 of file rgamma.c.

34 {
35 /* Check parameters */
36 if(x==0.0) return(0.0);
37 if((x<0.0) || (a<=0.0)) return(nan(""));
38
39 if((x>1.0) && (x>a)) return(1.0 - igamc(a, x));
40
41 /* Left tail of incomplete Gamma function:
42 x^a * e^-x * Sum(k=0,Inf)(x^k / Gamma(a+k+1)) */
43
44 double ans, ax, c, r;
45
46 /* Compute x**a * exp(-x) / Gamma(a) */
47 ax=a*log(x) - x - lgamma(a);
48 if(ax<-DBL_MAX_10_EXP) return(0.0); // underflow
49 ax=exp(ax);
50
51 /* power series */
52 r=a; c=1.0; ans=1.0;
53 do {
54 r+=1.0;
55 c*=x/r;
56 ans+=c;
57 } while(c/ans > DBL_EPSILON);
58 return(ans*ax/a);
59}
igamc
double igamc(double a, double x)
Definition rgamma.c:70

Referenced by igamc(), and mfEvalY().

◆ igamc()

double igamc ( double a,
double x )
extern

Regularized gamma function, more specifically, upper incomplete gamma function divided by gamma function.

f(a,x) = (1/Gamma(a)) * Integral(x,Inf)(e^-t * t^(a-1))dt Standard gamma distribution is assumed (Beta=1).

Returns
Returns the value of regularized gamma function, or NaN in case of an error.
See also
igam, inverfc, factorial
Parameters
aShape parameter alpha; must be > 0.
xIntegral from x to infinity; must be >= 0

Definition at line 70 of file rgamma.c.

75 {
76 /* Check parameters */
77 if((x<0.0) || (a<=0.0)) return(nan(""));
78
79 if((x<1.0) || (x<a)) return(1.0-igam(a, x));
80
81 double ans, ax, c, yc, r, t, y, z;
82 double pk, pkm1, pkm2, qk, qkm1, qkm2;
83 double big=4.503599627370496E+015;
84 double biginv=2.22044604925031308085E-016;
85
86 ax = a*log(x) - x - lgamma(a);
87 if(ax < -DBL_MAX_10_EXP) return(0.0); // underflow
88 ax=exp(ax);
89
90 /* continued fraction */
91 y=1.0-a; z=x+y+1.0; c=0.0;
92 pkm2=1.0; qkm2=x; pkm1=x+1.0; qkm1=z*x;
93 ans=pkm1/qkm1;
94 do {
95 c+=1.0; y+=1.0; z+=2.0;
96 yc=y*c; pk=pkm1*z - pkm2*yc; qk=qkm1*z - qkm2*yc;
97 if(qk!=0.0) {r=pk/qk; t=fabs((ans-r)/r); ans=r;}
98 else t=1.0;
99 pkm2=pkm1; pkm1=pk; qkm2=qkm1; qkm1=qk;
100 if(fabs(pk)>big) {pkm2*=biginv; pkm1*=biginv; qkm2*=biginv; qkm1*=biginv;}
101 } while(t>DBL_EPSILON);
102 return(ans*ax);
103}
igam
double igam(double a, double x)
Definition rgamma.c:29

Referenced by igam().

◆ mfEval2ndInt()

int mfEval2ndInt ( const char * fid,
const int parNr,
const double * p,
const int sampleNr,
const double * x,
double * v,
const int verbose )
extern

Evaluate the 2nd integral of function from 0 to specified x values.

Returns
Returns non-zero in case of an error.
Todo
Add and test more functions.
Author
Vesa Oikonen
See also
mfEvalInt, mfEvalY, mfEvalFrameY
Parameters
fidFunction code.
parNrNr of function parameters.
pArray of function parameters.
sampleNrSize of x and y arrays.
xArray of x values where function values are to be evaluated.
vPointer to allocated array of 2nd integral values where results will be written.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 700 of file func.c.

715 {
716 if(verbose>0) {
717 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
718 fflush(stdout);
719 }
720 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || v==NULL) return(1);
721
722 int i;
723 double xt, delayt=0.0;
724
725 if(strcasecmp(fid, "fengm2")==0) {
726 if(parNr<6) return(2);
727 if(parNr>6) delayt=p[6];
728 double A1, A2, A3, L1, L2, L3;
729 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
730 if(fabs(L1)<1.0E-06 || fabs(L2)<1.0E-06 || fabs(L3)<1.0E-06) {
731 if(verbose>1) printf("too small L parameter(s)\n");
732 return(3);
733 }
734 for(i=0; i<sampleNr; i++) {
735 xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
736 if(L1!=0.0) {
737 v[i] += (((A1*xt-A2-A3)*L1 - 2.0*A1) * exp(L1*xt) + 2.0*A1)/(L1*L1*L1);
738 v[i] += (A1*xt+A2+A3)/(L1*L1);
739 v[i] += (A2*xt+A3*xt)/L1;
740 }
741 if(L2!=0.0) v[i] -= (A2/(L2*L2))*(1.0 - exp(L2*xt)) + A2*xt/L2;
742 if(L3!=0.0) v[i] -= (A3/(L3*L3))*(1.0 - exp(L3*xt)) + A3*xt/L3;
743 }
744 return(0);
745 }
746 if(strcasecmp(fid, "fengm2e")==0) {
747 if(parNr<8) return(2);
748 if(parNr>8) delayt=p[8];
749 double A1, A2, A3, A4, L1, L2, L3, L4;
750 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
751 if(fabs(L1)<1.0E-15 || fabs(L2)<1.0E-15 || fabs(L3)<1.0E-15 || fabs(L4)<1.0E-15) {
752 if(verbose>1) printf("too small L parameter(s)\n");
753 return(3);
754 }
755 for(i=0; i<sampleNr; i++) {
756 xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
757 v[i] = A4*L1*L1*L1*L2*L2*L3*L3*L4*xt +
758 A4*L1*L1*L1*L2*L2*L3*L3*(1-exp(L4*xt)) +
759 L4*L4*(A3*L1*L1*L1*L2*L2*L3*xt + A3*L1*L1*L1*L2*L2*(1.0-exp(L3*xt)) +
760 L3*L3*(A2*L1*L1*L1*L2*xt + A2*L1*L1*L1*(1.0-exp(L2*xt)) -
761 L2*L2*((A2+A3+A4)*xt*L1*L1 + (A1*xt+A2+A3+A4)*L1 +
762 ((A1*xt-A2-A3-A4)*L1-2.0*A1)*exp(L1*xt) + 2.0*A1)));
763 v[i] /= -(L1*L1*L1*L2*L2*L3*L3*L4*L4);
764 }
765 return(0);
766 }
767 /* Gamma variate function when p[1]==1 */
768 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
769 if(parNr<3) return(2);
770 if(parNr>3) delayt=p[3];
771 for(i=0; i<sampleNr; i++) {
772 xt=x[i]-delayt; v[i]=0.0;
773 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
774 v[i] = p[0]*p[2]*p[2]*((p[2]+p[2]+xt)*exp(-xt/p[2]) - p[2] - p[2] + xt);
775 }
776 return(0);
777 }
778
779 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
780 return(10);
781}

◆ mfEvalFrameY()

int mfEvalFrameY ( const char * fid,
const int parNr,
const double * p,
const int sampleNr,
const double * x1,
const double * x2,
double * y,
const int verbose )
extern

Evaluate the PET time frame mean of function values.

Returns
Returns non-zero in case of an error.
Author
Vesa Oikonen
Todo
Only MF_FENGM2, MF_FENGM2E, and MF_DMSURGE are functional and tested.
Parameters
fidFunction code.
parNrNr of function parameters.
pArray of function parameters.
sampleNrSize of x and y arrays.
x1Array of frame start times.
x2Array of frame end times.
yPointer to allocated array of y values where function values will be written.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 790 of file func.c.

807 {
808 if(verbose>0) {
809 printf("%s('%s', %d, p[], %d, x1[], x2[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
810 fflush(stdout);
811 }
812 if(parNr<1 || p==NULL || sampleNr<1 || x1==NULL || x2==NULL || y==NULL) return(1);
813 if(verbose>10) {
814 printf("p[] := %g", p[0]);
815 for(int i=1; i<parNr; i++) printf(", %g", p[i]);
816 printf("\n");
817 }
818
819 int i, ret;
820 double xt1, xt2, v, fd, delayt=0.0;
821
822 if(strcasecmp(fid, "fengm2")==0) {
823 if(parNr<6) return(2);
824 if(parNr>6) delayt=p[6];
825 double A1, A2, A3, L1, L2, L3;
826 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
827 if(fabs(L1)<1.0E-12 || fabs(L2)<1.0E-12 || fabs(L3)<1.0E-12) {
828 if(verbose>1) printf("too small L parameter(s)\n");
829 return(3);
830 }
831 for(i=0; i<sampleNr; i++) {
832 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
833 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
834 if(!(xt2>0.0)) continue;
835 fd=xt2-xt1; // requested frame duration
836 if(fd<1.0E-025) {
837 // if frame duration is about zero, then use the method for that
838 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
839 if(ret!=0) return(ret);
840 continue;
841 }
842 if(!(xt1>0.0)) xt1=0.0; // integral start time must not be <0
843 /* Calculate integral between xt1 and xt2 */
844 y[i] = A3*L1*L1*L2*(exp(L3*xt1) - exp(L3*xt2)) +
845 L3*( A2*L1*L1*(exp(L2*xt1) - exp(L2*xt2)) +
846 L2*(((A1*xt1-A2-A3)*L1-A1)*exp(L1*xt1) -
847 ((A1*xt2-A2-A3)*L1-A1)*exp(L1*xt2)
848 )
849 );
850 y[i] /= -(L1*L1*L2*L3);
851 /* Divide integral with the original frame duration */
852 y[i]/=fd;
853 }
854 return(0);
855 }
856
857 if(strcasecmp(fid, "fengm2e")==0) {
858 if(parNr<8) return(2);
859 if(parNr>8) delayt=p[8];
860 double A1, A2, A3, A4, L1, L2, L3, L4;
861 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
862 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20 || fabs(L4)<1.0E-20) {
863 if(verbose>1) printf("too small L parameter(s)\n");
864 return(3);
865 }
866 for(i=0; i<sampleNr; i++) {
867 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
868 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
869 if(!(xt2>0.0)) continue;
870 fd=xt2-xt1; // requested frame duration
871 if(fd<1.0E-025) {
872 // if frame duration is about zero, then use the method for that
873 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
874 if(ret!=0) return(ret);
875 continue;
876 }
877 if(!(xt1>0.0)) xt1=0.0; // integral start time must not be <0
878 /* Calculate integral between xt1 and xt2 */
879 y[i] = A4*L1*L1*L2*L3*(exp(L4*xt1) - exp(L4*xt2)) +
880 L4*( A3*L1*L1*L2*(exp(L3*xt1) - exp(L3*xt2)) +
881 L3*( A2*L1*L1*(exp(L2*xt1) - exp(L2*xt2)) +
882 L2*(((A1*xt1-A2-A3-A4)*L1-A1)*exp(L1*xt1) -
883 ((A1*xt2-A2-A3-A4)*L1-A1)*exp(L1*xt2)
884 )));
885 y[i] /= -(L1*L1*L2*L3*L4);
886 /* Divide integral with the original frame duration */
887 y[i]/=fd;
888 }
889 return(0);
890 }
891 /* Gamma variate function when p[1]==1 */
892 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
893 if(parNr<3) return(2);
894 if(parNr>3) delayt=p[3];
895 for(i=0; i<sampleNr; i++) {
896 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
897 if(fabs(p[2])<1.0E-10) continue;
898 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
899 if(!(xt2>0.0)) continue;
900 fd=xt2-xt1; // requested frame duration
901 if(fd<1.0E-025) {
902 // if frame duration is about zero, then use the method for that
903 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
904 if(ret!=0) return(ret);
905 continue;
906 }
907 if(!(xt1>=0.0)) xt1=0.0; // integral start time must not be <0
908 /* Calculate integral between xt1 and xt2 */
909 y[i] = p[0]*p[2]*( (p[2]+xt1)*exp(-xt1/p[2]) - (p[2]+xt2)*exp(-xt2/p[2]) );
910 /* Divide integral with the original frame duration */
911 y[i]/=fd;
912 }
913 return(0);
914 }
915
916
917 /* Sum of Surge functions with delay, in form
918 f(x)=p1*(x-p0)*exp(-p2*(x-p0)) + p3*(x-p0)*exp(-p4*(x-p0)) + ... */
919 if(strcasecmp(fid, "dmsurge")==0) {
920 if(parNr<3) return(2);
921 delayt=p[0];
922 for(int i=0; i<sampleNr; i++) {
923 y[i]=0.0;
924 fd=x2[i]-x1[i]; // requested frame duration
925 if(verbose>11) printf(" fd[%d]=%g\n", i, fd);
926 if(!(fd>=0.0)) continue;
927 xt1=x1[i]-delayt; xt2=x2[i]-delayt;
928 if(!(xt2>0.0)) continue;
929 if(xt1<0.0) xt1=0.0; // integral start time must not be <0
930 if(fd<1.0E-025) {
931 // if frame duration is about zero, then simply calculate f(x) for that
932 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
933 if(ret!=0) return(ret);
934 if(verbose>11) printf(" y[%d]=%g\n", i, y[i]);
935 continue;
936 }
937 /* Calculate integral between xt1 and xt2 */
938 for(int pi=2; pi<parNr; pi+=2) {
939 double a, b;
940 a=p[pi-1]; b=p[pi];
941 if(!(a>0.0)) {
942 } else if(fabs(b)<1.0E-10) { // exp(0)=1
943 y[i] += a*(xt2*xt2 - xt1*xt1);
944 } else { // divider (b*b) must not be zero
945 y[i] += a*((b*xt1+1.0)*exp(-b*xt1) - (b*xt2+1.0)*exp(-b*xt2))/(b*b);
946 }
947 if(verbose>11) printf(" a=%g b=%g xt1=%g xt2=%g y[%d]=%g\n", a, b, xt1, xt2, i, y[i]);
948 }
949 /* Divide integral with the original frame duration */
950 y[i]/=fd;
951 if(verbose>11) printf(" y[%d]=%g\n", i, y[i]);
952 }
953 return(0);
954 }
955
956
957 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
958 return(10);
959}
mfEvalY
int mfEvalY(const char *fid, const int parNr, const double *p, const int sampleNr, const double *x, double *y, const int verbose)
Definition func.c:26

Referenced by spectralDMSurge().

◆ mfEvalInt()

int mfEvalInt ( const char * fid,
const int parNr,
const double * p,
const int sampleNr,
const double * x,
double * v,
const int verbose )
extern

Evaluate the function integral from 0 to specified x values.

Returns
Returns non-zero in case of an error.
Author
Vesa Oikonen
Todo
Add more functions and tests.
See also
mfEvalY, mfEvalFrameY, mfEval2ndInt, mfEvalIntToInf
Parameters
fidFunction code.
parNrNr of function parameters.
pArray of function parameters.
sampleNrSize of x and y arrays.
xArray of x values where function values are to be evaluated.
vPointer to allocated array of integral values where results will be written.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 436 of file func.c.

451 {
452 if(verbose>0) {
453 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
454 fflush(stdout);
455 }
456 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || v==NULL) return(1);
457
458 if(strcasecmp(fid, "fengm2")==0) {
459 if(parNr<6) return(2);
460 double delayt=0.0; if(parNr>6) delayt=p[6];
461 double a;
462 double A1, A2, A3, L1, L2, L3; // lambdas are negative
463 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
464 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20) {
465 if(verbose>1) printf("too small L parameter(s)\n");
466 return(3);
467 }
468 for(int i=0; i<sampleNr; i++) {
469 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
470 if(L1!=0.0) {
471 a=exp(L1*xt);
472 v[i] += (A1+L1*(A2+A3))*(1.0-a)/(L1*L1);
473 v[i] += (A1/L1)*xt*a;
474 }
475 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
476 if(L3!=0.0) v[i]-=(A3/L3)*(1.0-exp(L3*xt)); else v[i]+=A3*xt;
477 }
478 return(0);
479 }
480 if(strcasecmp(fid, "fengm2e")==0) {
481 if(parNr<8) return(2);
482 double delayt=0.0; if(parNr>8) delayt=p[8];
483 double a;
484 double A1, A2, A3, A4, L1, L2, L3, L4; // lambdas are negative
485 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
486 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20 || fabs(L4)<1.0E-20) {
487 if(verbose>1) printf("too small L parameter(s)\n");
488 return(3);
489 }
490 for(int i=0; i<sampleNr; i++) {
491 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
492 if(L1!=0.0) {
493 a=exp(L1*xt);
494 v[i] += (A1+L1*(A2+A3+A4))*(1.0-a)/(L1*L1);
495 v[i] += (A1/L1)*xt*a;
496 }
497 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
498 if(L3!=0.0) v[i]-=(A3/L3)*(1.0-exp(L3*xt)); else v[i]+=A3*xt;
499 if(L4!=0.0) v[i]-=(A4/L4)*(1.0-exp(L4*xt)); else v[i]+=A4*xt;
500 }
501 return(0);
502 }
503 if(strcasecmp(fid, "fengm2s")==0) {
504 if(parNr<4) return(2);
505 double delayt=0.0; if(parNr>4) delayt=p[4];
506 double a;
507 double A1, A2, L1, L2; // lambdas are negative
508 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3];
509 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20) {
510 if(verbose>1) printf("too small L parameter(s)\n");
511 return(3);
512 }
513 for(int i=0; i<sampleNr; i++) {
514 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
515 if(L1!=0.0) {
516 a=exp(L1*xt);
517 v[i] += (A1+L1*A2)*(1.0-a)/(L1*L1);
518 v[i] += (A1/L1)*xt*a;
519 }
520 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
521 }
522 return(0);
523 }
524 /* Gamma variate function when p[1]==1 */
525 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
526 if(parNr<3) return(2);
527 double delayt=0.0; if(parNr>3) delayt=p[3];
528 for(int i=0; i<sampleNr; i++) {
529 double xt=x[i]-delayt; v[i]=0.0;
530 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
531 v[i] = p[0]*p[2]*(p[2] - (p[2]+xt)*exp(-xt/p[2]));
532 }
533 return(0);
534 }
535
536 /* Integral of sum of Surge functions with delay, in form
537 f(x)=p1*(x-p0)*exp(-p2*(x-p0)) + p3*(x-p0)*exp(-p4*(x-p0)) + ... */
538 if(strcasecmp(fid, "dmsurge")==0) {
539 if(parNr<3) return(2);
540 for(int i=0; i<sampleNr; i++) {
541 double xt=x[i]-p[0];
542 v[i]=0.0;
543 if(xt>0.0)
544 for(int pi=2; pi<parNr; pi+=2) {
545 double a, b;
546 a=p[pi-1]; b=p[pi];
547 if(a>0.0) {
548 if(fabs(b)>1.0E-12)
549 v[i] += a*(1.0-(b*xt+1.0)*exp(-b*xt))/(b*b);
550 else
551 v[i] += a*xt*xt;
552 }
553 }
554 }
555 return(0);
556 }
557
558 /* Surge function in form f(x)=a*x*exp(-b*x)*b^2, to give AUC=a from 0 to infinity. */
559 if(strcasecmp(fid, "surge")==0) {
560 if(parNr<2) return(2);
561 for(int i=0; i<sampleNr; i++) {
562 double xt=x[i]; v[i]=0.0;
563 if(!(xt>0.0)) continue;
564 v[i] = p[0]*(1.0 - (p[1]*xt + 1.0)*exp(-p[1]*xt));
565 }
566 return(0);
567 }
568
569 /* Surge function in form f(x)=a*x*exp(-b*x). Maximum value is at 1/b. */
570 if(strcasecmp(fid, "tradsurge")==0) {
571 if(parNr<2) return(2);
572 double f=p[0]/(p[1]*p[1]);
573 for(int i=0; i<sampleNr; i++) {
574 double xt=x[i]; v[i]=0.0;
575 if(!(xt>0.0)) continue;
576 v[i] = f*(1.0 - (p[1]*xt + 1.0)*exp(-p[1]*xt));
577 }
578 return(0);
579 }
580
581 /* Surge function for late FDG AIF with time delay */
582 if(strcasecmp(fid, "surgefdgaif")==0) {
583 if(parNr<5) return(2);
584 double delayt=p[0];
585 for(int i=0; i<sampleNr; i++) {
586 double xt=x[i]-delayt; v[i]=0.0;
587 if(!(xt>0.0)) continue;
588 v[i] += p[3]*(1.0-(p[4]*p[1]*xt + 1.0)*exp(-p[4]*p[1]*xt))/(p[1]*p[1]*p[4]*p[4]);
589 v[i] += (1.0-exp(-p[1]*xt))/p[1];
590 v[i] *= p[2];
591 }
592 return(0);
593 }
594
595 /* Probability density function of Erlang distribution */
596 if(strcasecmp(fid, "erlangpdf")==0 && p[2]<20.5) { // Supports only N=0,1,2,...,20
597 if(parNr<3) return(2);
598 double A=p[0], k=p[1]; if(!(k>=0.0)) return(2);
599 int N=(int)round(p[2]); if(!(N>=0) || N>20) return(2);
600 if(N==0)
601 for(int i=0; i<sampleNr; i++) {
602 if(x[i]>=0.0) v[i]=A; else v[i]=0.0;
603 }
604 else if(N==1)
605 for(int i=0; i<sampleNr; i++) {
606 if(x[i]>=0.0) v[i]=A*(1.0-exp(-k*x[i])); else v[i]=0.0;
607 }
608 else if(N==2)
609 for(int i=0; i<sampleNr; i++) {
610 if(x[i]>=0.0) v[i]=A*(1.0 - exp(-k*x[i]) - k*x[i]*exp(-k*x[i])); else v[i]=0.0;
611 }
612 else {
613 for(int i=0; i<sampleNr; i++) {
614 if(!(x[i]>=0.0)) {v[i]=0.0; continue;}
615 double s=1.0+k*x[i];
616 for(int j=2; j<N; j++) s+=pow(k*x[i], j)/lfactorial(j);
617 v[i]=A*(1.0-s*exp(-k*x[i]));
618 }
619 }
620
621 return(0);
622 }
623
624 /* Exponential functions */
625 if(strcasecmp(fid, "1exp")==0) {
626 if(parNr<2) return(2);
627 double A=p[0], k=p[1];
628 if(fabs(k)<1.0E-20) {
629 for(int i=0; i<sampleNr; i++) v[i]=A*x[i];
630 } else {
631 for(int i=0; i<sampleNr; i++) v[i]=(A/k)*(exp(k*x[i]) - 1.0);
632 }
633 return(0);
634 }
635 if(strcasecmp(fid, "2exp")==0) {
636 if(parNr<4) return(2);
637 for(int i=0; i<sampleNr; i++) v[i]=0.0;
638 for(int n=0; n<=1; n++) {
639 double A=p[n*2], k=p[n*2+1];
640 if(fabs(k)<1.0E-20) {
641 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
642 } else {
643 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
644 }
645 }
646 return(0);
647 }
648 if(strcasecmp(fid, "3exp")==0) {
649 if(parNr<6) return(2);
650 for(int i=0; i<sampleNr; i++) v[i]=0.0;
651 for(int n=0; n<=2; n++) {
652 double A=p[n*2], k=p[n*2+1];
653 if(fabs(k)<1.0E-20) {
654 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
655 } else {
656 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
657 }
658 }
659 return(0);
660 }
661 if(strcasecmp(fid, "4exp")==0) {
662 if(parNr<8) return(2);
663 for(int i=0; i<sampleNr; i++) v[i]=0.0;
664 for(int n=0; n<=3; n++) {
665 double A=p[n*2], k=p[n*2+1];
666 if(fabs(k)<1.0E-20) {
667 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
668 } else {
669 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
670 }
671 }
672 return(0);
673 }
674 if(strcasecmp(fid, "5exp")==0) {
675 if(parNr<10) return(2);
676 for(int i=0; i<sampleNr; i++) v[i]=0.0;
677 for(int n=0; n<=4; n++) {
678 double A=p[n*2], k=p[n*2+1];
679 if(fabs(k)<1.0E-20) {
680 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
681 } else {
682 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
683 }
684 }
685 return(0);
686 }
687
688 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
689 return(10);
690}
lfactorial
unsigned long long int lfactorial(unsigned long long int n)
Definition intutil.c:63

◆ mfEvalIntToInf()

int mfEvalIntToInf ( const char * fid,
const int parNr,
const double * p,
const double x,
double * v,
const int verbose )
extern

Evaluate the function integral from specified x to infinity.

Returns
Returns non-zero in case of an error.
Author
Vesa Oikonen
Todo
Only exponential functions are functional and tested.
See also
mfEvalInt, mfEvalY
Parameters
fidFunction code.
parNrNr of function parameters.
pArray of function parameters.
xX value; must be 0 or larger.
vPointer for result integral.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 969 of file func.c.

982 {
983 if(verbose>0) {
984 printf("%s('%s', %d, p[], %g, y, %d)\n", __func__, fid, parNr, x, verbose);
985 fflush(stdout);
986 }
987 if(parNr<1 || p==NULL || v==NULL) return(1);
988 if(!(x>=0.0)) return(2);
989 *v=0.0;
990
991 /* Exponential functions */
992 if(strcasecmp(fid, "1exp")==0) {
993 if(parNr<2) return(2);
994 double A=p[0], k=p[1];
995 *v=(-A/k)*exp(k*x);
996 return(0);
997 }
998 if(strcasecmp(fid, "2exp")==0) {
999 if(parNr<4) return(2);
1000 for(int n=0; n<=1; n++) {
1001 double A=p[n*2], k=p[n*2+1];
1002 *v+=(-A/k)*exp(k*x);
1003 }
1004 return(0);
1005 }
1006 if(strcasecmp(fid, "3exp")==0) {
1007 if(parNr<6) return(2);
1008 for(int n=0; n<=2; n++) {
1009 double A=p[n*2], k=p[n*2+1];
1010 *v+=(-A/k)*exp(k*x);
1011 }
1012 return(0);
1013 }
1014 if(strcasecmp(fid, "4exp")==0) {
1015 if(parNr<8) return(2);
1016 for(int n=0; n<=3; n++) {
1017 double A=p[n*2], k=p[n*2+1];
1018 *v+=(-A/k)*exp(k*x);
1019 }
1020 return(0);
1021 }
1022 if(strcasecmp(fid, "5exp")==0) {
1023 if(parNr<10) return(2);
1024 for(int n=0; n<=4; n++) {
1025 double A=p[n*2], k=p[n*2+1];
1026 *v+=(-A/k)*exp(k*x);
1027 }
1028 return(0);
1029 }
1030
1031
1032
1033 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
1034 return(10);
1035}

◆ mfEvalY()

int mfEvalY ( const char * fid,
const int parNr,
const double * p,
const int sampleNr,
const double * x,
double * y,
const int verbose )
extern

Evaluate the function values at specified x values.

Returns
Returns non-zero in case of an error.
Author
Vesa Oikonen
Todo
Add support for more functions.
See also
mfEvalInt, mfEvalFrameY, parExampleTTACs, parExamplePerfectBolus, mfCreateTAC
Parameters
fidFunction code.
parNrNr of function parameters.
pArray of function parameters.
sampleNrSize of x and y arrays.
xArray of x values where function values are to be evaluated.
yPointer to allocated array of y values where function values will be written.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 26 of file func.c.

41 {
42 if(verbose>0) {
43 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
44 fflush(stdout);
45 }
46 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || y==NULL) return(1);
47
48 if(strcasecmp(fid, "step")==0) {
49 /* Parameter number must be an even number and >=2 */
50 if(parNr<2 || parNr&1) return(2);
51 for(int i=0; i<sampleNr; i++) {
52 y[i]=0.0;
53 for(int j=0; j<parNr; j+=2) if(x[i]>=p[j]) y[i]=p[j+1];
54 }
55 return(0);
56 }
57
58 if(strcasecmp(fid, "level")==0) {
59 for(int i=0; i<sampleNr; i++) y[i]=p[0];
60 return(0);
61 }
62
63 if(strcasecmp(fid, "fengm2")==0) {
64 if(parNr<6) return(2);
65 double delayt=0.0; if(parNr>6) delayt=p[6];
66 for(int i=0; i<sampleNr; i++) {
67 double xt=x[i]-delayt; y[i]=0.0; if(!(xt>0.0)) continue;
68 y[i] += (p[0]*xt-p[2]-p[4])*exp(p[1]*xt);
69 y[i] += p[2]*exp(p[3]*xt);
70 y[i] += p[4]*exp(p[5]*xt);
71 }
72 return(0);
73 }
74 if(strcasecmp(fid, "fengm2s")==0) {
75 if(parNr<4) return(2);
76 double delayt=0.0; if(parNr>4) delayt=p[4];
77 for(int i=0; i<sampleNr; i++) {
78 double xt=x[i]-delayt; y[i]=0.0; if(!(xt>0.0)) continue;
79 y[i] += (p[0]*xt-p[2])*exp(p[1]*xt);
80 y[i] += p[2]*exp(p[3]*xt);
81 }
82 return(0);
83 }
84 if(strcasecmp(fid, "fengm2e")==0) {
85 if(parNr<8) return(2);
86 double delayt=0.0; if(parNr>8) delayt=p[8];
87 for(int i=0; i<sampleNr; i++) {
88 double xt=x[i]-delayt; y[i]=0.0; if(!(xt>0.0)) continue;
89 y[i] += (p[0]*xt-p[2]-p[4]-p[6])*exp(p[1]*xt);
90 y[i] += p[2]*exp(p[3]*xt);
91 y[i] += p[4]*exp(p[5]*xt);
92 y[i] += p[6]*exp(p[7]*xt);
93 }
94 return(0);
95 }
96 /* Gamma variate function */
97 if(strcasecmp(fid, "gammav")==0) {
98 if(parNr<3) return(2);
99 double delayt=0.0; if(parNr>3) delayt=p[3];
100 for(int i=0; i<sampleNr; i++) {
101 double xt=x[i]-delayt; y[i]=0.0;
102 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
103 y[i] = p[0]*pow(xt, p[1])*exp(-xt/p[2]);
104 }
105 return(0);
106 }
107 /* Exponential function */
108 if(strcasecmp(fid, "1exp")==0) {
109 if(parNr<2) return(2);
110 for(int i=0; i<sampleNr; i++) {
111 y[i]=0.0;
112 y[i] += p[0]*exp(p[1]*x[i]);
113 }
114 return(0);
115 }
116 /* Sum of two exponentials */
117 if(strcasecmp(fid, "2exp")==0) {
118 if(parNr<4) return(2);
119 for(int i=0; i<sampleNr; i++) {
120 y[i]=0.0;
121 y[i] += p[0]*exp(p[1]*x[i]);
122 y[i] += p[2]*exp(p[3]*x[i]);
123 }
124 return(0);
125 }
126 /* Sum of three exponentials */
127 if(strcasecmp(fid, "3exp")==0) {
128 if(parNr<6) return(2);
129 for(int i=0; i<sampleNr; i++) {
130 y[i]=0.0;
131 y[i] += p[0]*exp(p[1]*x[i]);
132 y[i] += p[2]*exp(p[3]*x[i]);
133 y[i] += p[4]*exp(p[5]*x[i]);
134 }
135 return(0);
136 }
137 /* Sum of four exponentials */
138 if(strcasecmp(fid, "4exp")==0) {
139 if(parNr<8) return(2);
140 for(int i=0; i<sampleNr; i++) {
141 y[i]=0.0;
142 y[i] += p[0]*exp(p[1]*x[i]);
143 y[i] += p[2]*exp(p[3]*x[i]);
144 y[i] += p[4]*exp(p[5]*x[i]);
145 y[i] += p[6]*exp(p[7]*x[i]);
146 }
147 return(0);
148 }
149 /* Sum of five exponentials */
150 if(strcasecmp(fid, "5exp")==0) {
151 if(parNr<10) return(2);
152 for(int i=0; i<sampleNr; i++) {
153 y[i]=0.0;
154 y[i] += p[0]*exp(p[1]*x[i]);
155 y[i] += p[2]*exp(p[3]*x[i]);
156 y[i] += p[4]*exp(p[5]*x[i]);
157 y[i] += p[6]*exp(p[7]*x[i]);
158 y[i] += p[8]*exp(p[9]*x[i]);
159 }
160 return(0);
161 }
162
163 /* Inverted gamma cdf for plasma parent fraction */
164 if(strcasecmp(fid, "ppfigamc")==0) {
165 if(parNr<4) return(2);
166 double delayt=0.0; if(parNr>4) delayt=p[4];
167 double a=p[0];
168 double b=p[1];
169 double c=p[2]; if(c<0.0) return(2);
170 double d=p[3]; if(d<=0.0) return(2);
171 for(int i=0; i<sampleNr; i++) {
172 double t=x[i]-delayt; if(t<0.0) t=0.0;
173 y[i]=a*(1.0-b*igam(d, c*t));
174 }
175 return(0);
176 }
177
178 /* Weibull cdf with time delay */
179 if(strcasecmp(fid, "weibullcdfd")==0) {
180 if(parNr<4) return(2);
181 double delayt=p[0];
182 for(int i=0; i<sampleNr; i++) {
183 double xt=x[i]-delayt; y[i]=0.0;
184 if(!(xt>0.0)) continue;
185 y[i] = p[1]*(1.0-exp(-pow(xt/p[2], p[3])));
186 }
187 return(0);
188 }
189
190 /* Weibull cdf with time delay and cdf derivative */
191 if(strcasecmp(fid, "weibullcdfdd")==0) {
192 if(parNr<5) return(2);
193 double delayt=p[0];
194 for(int i=0; i<sampleNr; i++) {
195 double xt=x[i]-delayt; y[i]=0.0;
196 if(!(xt>0.0)) continue;
197 double a=exp(-pow(xt/p[2], p[3]));
198 double b=(p[3]/p[2]) * pow(xt/p[2], p[3]-1.0) * a;
199 y[i] = p[1]*(b + p[4]*(1.0-a));
200 }
201 return(0);
202 }
203
204 /* Sum of Surge functions with delay, in form
205 f(x)=p1*(x-p0)*exp(-p2*(x-p0)) + p3*(x-p0)*exp(-p4*(x-p0)) + ... */
206 if(strcasecmp(fid, "dmsurge")==0) {
207 if(parNr<3) return(2);
208 double delayt=p[0];
209 for(int i=0; i<sampleNr; i++) {
210 double xt=x[i]-delayt;
211 y[i]=0.0;
212 if(xt>0.0)
213 for(int pi=2; pi<parNr; pi+=2) {
214 double a, b;
215 a=p[pi-1]; b=p[pi];
216 y[i] += a*xt*exp(-b*xt);
217 }
218 }
219 return(0);
220 }
221
222 /* Surge function in form f(x)=a*x*exp(-b*x)*b^2 to give AUC=a from 0 to infinity.
223 Maximum value is at 1/b. */
224 if(strcasecmp(fid, "surge")==0) {
225 if(parNr<2) return(2);
226 double f=p[0]*(p[1]*p[1]);
227 for(int i=0; i<sampleNr; i++) {
228 double xt=x[i]; y[i]=0.0;
229 if(!(xt>0.0)) continue;
230 y[i] = f*xt*exp(-p[1]*xt);
231 }
232 return(0);
233 }
234
235 /* Surge function in form f(x)=a*x*exp(-b*x). Maximum value is at 1/b. */
236 if(strcasecmp(fid, "tradsurge")==0) {
237 if(parNr<2) return(2);
238 for(int i=0; i<sampleNr; i++) {
239 double xt=x[i]; y[i]=0.0;
240 if(!(xt>0.0)) continue;
241 y[i] = p[0]*xt*exp(-p[1]*xt);
242 }
243 return(0);
244 }
245
246 /* Surge function with recirculation, in form
247 f(x)=a*(x*exp(-b*x) + (c/(b*b))*(1-(b*x+1)*exp(-b*x))). */
248 if(strcasecmp(fid, "surgerecirc")==0) {
249 if(parNr<2) return(2);
250 double c=0.0;
251 if(parNr>=3) c=p[2];
252 for(int i=0; i<sampleNr; i++) {
253 double xt=x[i]; y[i]=0.0;
254 if(!(xt>0.0)) continue;
255 double e=exp(-p[1]*xt);
256 y[i] = xt*e + (c/(p[1]*p[1]))*(1.0 - (p[1]*xt+1.0)*e);
257 y[i] *= p[0];
258 }
259 return(0);
260 }
261
262 /* Surge function with recirculation and delay, in form
263 f(x)=a*((x-d)*exp(-b*(x-d)) + (c/(b*b))*(1-(b*(x-d)+1)*exp(-b*(x-d)))). */
264 if(strcasecmp(fid, "surgerecircd")==0) {
265 if(parNr<2) return(2);
266 double c=0.0; if(parNr>=3) c=p[2];
267 double delayt=0.0; if(parNr>=4) delayt=p[3];
268 for(int i=0; i<sampleNr; i++) {
269 double xt=x[i]-delayt; y[i]=0.0;
270 if(!(xt>0.0)) continue;
271 double e=exp(-p[1]*xt);
272 y[i] = xt*e + (c/(p[1]*p[1]))*(1.0 - (p[1]*xt+1.0)*e);
273 y[i] *= p[0];
274 }
275 return(0);
276 }
277
278 /* Surge function with recirculation, for plasma-to-blood ratio.
279 RBC-to-plasma r(x)=a*(x*exp(-b*x) + (c/(b*b))*(1-(b*x+1)*exp(-b*x)))
280 Plasma-to-blood f(x)=1/(1-HCT*(1-r(x))). */
281 if(strcasecmp(fid, "p2bsrc")==0) {
282 if(parNr<3) return(2);
283 double c=0.0; if(parNr>=4) c=p[3];
284 double HCT=p[0];
285 double a=p[1];
286 double b=p[2];
287 for(int i=0; i<sampleNr; i++) {
288 double xt=x[i]; if(!(xt>0.0)) {y[i]=1.0/(1.0-HCT); continue;}
289 double e=exp(-b*xt);
290 double r=a*(xt*e + (c/(b*b))*(1.0 - (b*xt+1.0)*e));
291 y[i] = 1.0/(1.0-HCT*(1.0-r));
292 }
293 return(0);
294 }
295
296 /* Feng M2 function for plasma-to-blood ratio.
297 RBC-to-plasma ratio r(x) using Feng M2, and then Plasma-to-blood f(x)=1/(1-HCT*(1-r(x))). */
298 if(strcasecmp(fid, "p2bfm2")==0) {
299 if(parNr<3) return(2);
300 double HCT=p[0];
301 double A1=p[1];
302 double L1=p[2];
303 double A2=0.0; if(parNr>3) A2=p[3];
304 double L2=0.0; if(parNr>4) L2=p[4];
305 double A3=0.0; if(parNr>5) A3=p[5];
306 double L3=0.0; if(parNr>6) L3=p[6];
307 for(int i=0; i<sampleNr; i++) {
308 double xt=x[i]; if(!(xt>0.0)) {y[i]=1.0/(1.0-HCT); continue;}
309 double r=(A1*xt - A2 - A3)*exp(-L1*xt) + A2*exp(-L2*xt) + A3*exp(-L3*xt);
310 y[i] = 1.0/(1.0-HCT*(1.0-r));
311 }
312 return(0);
313 }
314
315 /* Surge function for late FDG AIF with time delay */
316 if(strcasecmp(fid, "surgefdgaif")==0) {
317 if(parNr<5) return(2);
318 double delayt=p[0];
319 for(int i=0; i<sampleNr; i++) {
320 double xt=x[i]-delayt; y[i]=0.0;
321 if(!(xt>0.0)) continue;
322 y[i] = p[2]*(p[3]*xt*exp(-p[4]*p[1]*xt) + exp(-p[1]*xt));
323 }
324 return(0);
325 }
326
327 /* Probability density function of Erlang distribution */
328 if(strcasecmp(fid, "erlangpdf")==0) {
329 if(parNr<3) return(2);
330 double A=p[0], k=p[1]; if(!(k>=0.0)) return(2);
331 int N=(int)round(p[2]); if(!(N>0)) return(2);
332 unsigned long long int f=lfactorial(N-1);
333 for(int i=0; i<sampleNr; i++) {
334 y[i]=0.0; if(!(x[i]>=0.0)) continue;
335 y[i] = A*pow(k, N)*pow(x[i], N-1)*exp(-k*x[i])/f;
336 }
337 return(0);
338 }
339
340 /* Rational functions */
341 if(strcasecmp(fid, "ratf11")==0) {
342 if(parNr<4) return(2);
343 for(int i=0; i<sampleNr; i++) {
344 double a=p[0]+p[2]*x[i];
345 double b=p[1]+p[3]*x[i];
346 y[i]=a/b;
347 }
348 return(0);
349 }
350 if(strcasecmp(fid, "ratf21")==0) {
351 if(parNr<5) return(2);
352 for(int i=0; i<sampleNr; i++) {
353 double a=p[0]+p[2]*x[i]+p[4]*x[i]*x[i];
354 double b=p[1]+p[3]*x[i];
355 y[i]=a/b;
356 }
357 return(0);
358 }
359 if(strcasecmp(fid, "ratf22")==0) {
360 if(parNr<6) return(2);
361 for(int i=0; i<sampleNr; i++) {
362 double sqrx=x[i]*x[i];
363 double a=p[0]+p[2]*x[i]+p[4]*sqrx;
364 double b=p[1]+p[3]*x[i]+p[5]*sqrx;
365 y[i]=a/b;
366 }
367 return(0);
368 }
369 if(strcasecmp(fid, "ratf32")==0) {
370 if(parNr<7) return(2);
371 for(int i=0; i<sampleNr; i++) {
372 double sqrx=x[i]*x[i];
373 double a=p[0]+p[2]*x[i]+p[4]*sqrx+p[6]*sqrx*x[i];
374 double b=p[1]+p[3]*x[i]+p[5]*sqrx;
375 y[i]=a/b;
376 }
377 return(0);
378 }
379 if(strcasecmp(fid, "ratf33")==0) {
380 if(parNr<8) return(2);
381 for(int i=0; i<sampleNr; i++) {
382 double sqrx=x[i]*x[i];
383 double cubx=sqrx*x[i];
384 double a=p[0]+p[2]*x[i]+p[4]*sqrx+p[6]*cubx;
385 double b=p[1]+p[3]*x[i]+p[5]*sqrx+p[7]*cubx;
386 y[i]=a/b;
387 }
388 return(0);
389 }
390 if(strcasecmp(fid, "p2brf")==0) {
391 if(parNr<7) return(2);
392 for(int i=0; i<sampleNr; i++) {
393 if(x[i]<=0.0) {
394 y[i]=1.0/(1.0-p[0]);
395 } else {
396 double sqrx=x[i]*x[i];
397 double cubx=sqrx*x[i];
398 double a=0.0+p[1]*x[i]+p[3]*sqrx+p[5]*cubx;
399 double b=1.0+p[2]*x[i]+p[4]*sqrx+p[6]*cubx;
400 y[i]=1.0/(1.0-p[0]*(1.0-(a/b)));
401 }
402 }
403 return(0);
404 }
405
406 /* Hill functions */
407 if(strcasecmp(fid, "mpfhill")==0) {
408 if(parNr<3) return(2);
409 for(int i=0; i<sampleNr; i++) {
410 y[i]=p[0]*pow(x[i], p[1]) / (pow(x[i], p[1]) + p[2]);
411 }
412 return(0);
413 }
414 if(strcasecmp(fid, "p2bhill")==0) {
415 if(parNr<4) return(2);
416 for(int i=0; i<sampleNr; i++) {
417 double cpr=p[1]*pow(x[i], p[2]) / (pow(x[i], p[2]) + p[3]);
418 y[i]=1.0/(1.0-p[0]*(1.0-cpr));
419 }
420 return(0);
421 }
422
423
424 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
425 return(10);
426}

Referenced by mfCreateTAC(), mfEvalFrameY(), and spectralDMSurge().

Generated on for TPCCLIB by doxygen 1.15.0