Doxygen Source Code Documentation
cdf_12.c File Reference
#include "cdflib.h"
Go to the source code of this file.
Functions | |
double | brcomp (double *a, double *b, double *x, double *y) |
Function Documentation
|
Definition at line 2 of file cdf_12.c. References a, algdiv(), alnrel(), bcorr(), betaln(), brcomp(), c, fifdmax1(), fifdmin1(), gam1(), gamln1(), i, rlog1(), v, x0, and y0. Referenced by bfrac(), and brcomp().
00008 { 00009 static double Const = .398942280401433e0; 00010 static double brcomp,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z; 00011 static int i,n; 00012 /* 00013 ----------------- 00014 CONST = 1/SQRT(2*PI) 00015 ----------------- 00016 */ 00017 static double T1,T2; 00018 /* 00019 .. 00020 .. Executable Statements .. 00021 */ 00022 brcomp = 0.0e0; 00023 if(*x == 0.0e0 || *y == 0.0e0) return brcomp; 00024 a0 = fifdmin1(*a,*b); 00025 if(a0 >= 8.0e0) goto S130; 00026 if(*x > 0.375e0) goto S10; 00027 lnx = log(*x); 00028 T1 = -*x; 00029 lny = alnrel(&T1); 00030 goto S30; 00031 S10: 00032 if(*y > 0.375e0) goto S20; 00033 T2 = -*y; 00034 lnx = alnrel(&T2); 00035 lny = log(*y); 00036 goto S30; 00037 S20: 00038 lnx = log(*x); 00039 lny = log(*y); 00040 S30: 00041 z = *a*lnx+*b*lny; 00042 if(a0 < 1.0e0) goto S40; 00043 z -= betaln(a,b); 00044 brcomp = exp(z); 00045 return brcomp; 00046 S40: 00047 /* 00048 ----------------------------------------------------------------------- 00049 PROCEDURE FOR A .LT. 1 OR B .LT. 1 00050 ----------------------------------------------------------------------- 00051 */ 00052 b0 = fifdmax1(*a,*b); 00053 if(b0 >= 8.0e0) goto S120; 00054 if(b0 > 1.0e0) goto S70; 00055 /* 00056 ALGORITHM FOR B0 .LE. 1 00057 */ 00058 brcomp = exp(z); 00059 if(brcomp == 0.0e0) return brcomp; 00060 apb = *a+*b; 00061 if(apb > 1.0e0) goto S50; 00062 z = 1.0e0+gam1(&apb); 00063 goto S60; 00064 S50: 00065 u = *a+*b-1.e0; 00066 z = (1.0e0+gam1(&u))/apb; 00067 S60: 00068 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z; 00069 brcomp = brcomp*(a0*c)/(1.0e0+a0/b0); 00070 return brcomp; 00071 S70: 00072 /* 00073 ALGORITHM FOR 1 .LT. B0 .LT. 8 00074 */ 00075 u = gamln1(&a0); 00076 n = b0-1.0e0; 00077 if(n < 1) goto S90; 00078 c = 1.0e0; 00079 for(i=1; i<=n; i++) { 00080 b0 -= 1.0e0; 00081 c *= (b0/(a0+b0)); 00082 } 00083 u = log(c)+u; 00084 S90: 00085 z -= u; 00086 b0 -= 1.0e0; 00087 apb = a0+b0; 00088 if(apb > 1.0e0) goto S100; 00089 t = 1.0e0+gam1(&apb); 00090 goto S110; 00091 S100: 00092 u = a0+b0-1.e0; 00093 t = (1.0e0+gam1(&u))/apb; 00094 S110: 00095 brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t; 00096 return brcomp; 00097 S120: 00098 /* 00099 ALGORITHM FOR B0 .GE. 8 00100 */ 00101 u = gamln1(&a0)+algdiv(&a0,&b0); 00102 brcomp = a0*exp(z-u); 00103 return brcomp; 00104 S130: 00105 /* 00106 ----------------------------------------------------------------------- 00107 PROCEDURE FOR A .GE. 8 AND B .GE. 8 00108 ----------------------------------------------------------------------- 00109 */ 00110 if(*a > *b) goto S140; 00111 h = *a/ *b; 00112 x0 = h/(1.0e0+h); 00113 y0 = 1.0e0/(1.0e0+h); 00114 lambda = *a-(*a+*b)**x; 00115 goto S150; 00116 S140: 00117 h = *b/ *a; 00118 x0 = 1.0e0/(1.0e0+h); 00119 y0 = h/(1.0e0+h); 00120 lambda = (*a+*b)**y-*b; 00121 S150: 00122 e = -(lambda/ *a); 00123 if(fabs(e) > 0.6e0) goto S160; 00124 u = rlog1(&e); 00125 goto S170; 00126 S160: 00127 u = e-log(*x/x0); 00128 S170: 00129 e = lambda/ *b; 00130 if(fabs(e) > 0.6e0) goto S180; 00131 v = rlog1(&e); 00132 goto S190; 00133 S180: 00134 v = e-log(*y/y0); 00135 S190: 00136 z = exp(-(*a*u+*b*v)); 00137 brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b)); 00138 return brcomp; 00139 } /* END */ |