Doxygen Source Code Documentation
cdf_04.c File Reference
#include "cdflib.h"Go to the source code of this file.
Functions | |
| double | basym (double *a, double *b, double *lambda, double *eps) |
Function Documentation
|
||||||||||||||||||||
|
Definition at line 2 of file cdf_04.c. References a, basym(), bcorr(), c, erfc1(), i, r, rlog1(), z0, and zn. Referenced by basym(), and bratio().
00011 {
00012 static double e0 = 1.12837916709551e0;
00013 static double e1 = .353553390593274e0;
00014 static int num = 20;
00015 /*
00016 ------------------------
00017 ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
00018 ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
00019 THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
00020 ------------------------
00021 E0 = 2/SQRT(PI)
00022 E1 = 2**(-3/2)
00023 ------------------------
00024 */
00025 static int K3 = 1;
00026 static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
00027 z2,zn,znm1;
00028 static int i,im1,imj,j,m,mm1,mmj,n,np1;
00029 static double a0[21],b0[21],c[21],d[21],T1,T2;
00030 /*
00031 ..
00032 .. Executable Statements ..
00033 */
00034 basym = 0.0e0;
00035 if(*a >= *b) goto S10;
00036 h = *a/ *b;
00037 r0 = 1.0e0/(1.0e0+h);
00038 r1 = (*b-*a)/ *b;
00039 w0 = 1.0e0/sqrt(*a*(1.0e0+h));
00040 goto S20;
00041 S10:
00042 h = *b/ *a;
00043 r0 = 1.0e0/(1.0e0+h);
00044 r1 = (*b-*a)/ *a;
00045 w0 = 1.0e0/sqrt(*b*(1.0e0+h));
00046 S20:
00047 T1 = -(*lambda/ *a);
00048 T2 = *lambda/ *b;
00049 f = *a*rlog1(&T1)+*b*rlog1(&T2);
00050 t = exp(-f);
00051 if(t == 0.0e0) return basym;
00052 z0 = sqrt(f);
00053 z = 0.5e0*(z0/e1);
00054 z2 = f+f;
00055 a0[0] = 2.0e0/3.0e0*r1;
00056 c[0] = -(0.5e0*a0[0]);
00057 d[0] = -c[0];
00058 j0 = 0.5e0/e0*erfc1(&K3,&z0);
00059 j1 = e1;
00060 sum = j0+d[0]*w0*j1;
00061 s = 1.0e0;
00062 h2 = h*h;
00063 hn = 1.0e0;
00064 w = w0;
00065 znm1 = z;
00066 zn = z2;
00067 for(n=2; n<=num; n+=2) {
00068 hn = h2*hn;
00069 a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
00070 np1 = n+1;
00071 s += hn;
00072 a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
00073 for(i=n; i<=np1; i++) {
00074 r = -(0.5e0*((double)i+1.0e0));
00075 b0[0] = r*a0[0];
00076 for(m=2; m<=i; m++) {
00077 bsum = 0.0e0;
00078 mm1 = m-1;
00079 for(j=1; j<=mm1; j++) {
00080 mmj = m-j;
00081 bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
00082 }
00083 b0[m-1] = r*a0[m-1]+bsum/(double)m;
00084 }
00085 c[i-1] = b0[i-1]/((double)i+1.0e0);
00086 dsum = 0.0e0;
00087 im1 = i-1;
00088 for(j=1; j<=im1; j++) {
00089 imj = i-j;
00090 dsum += (d[imj-1]*c[j-1]);
00091 }
00092 d[i-1] = -(dsum+c[i-1]);
00093 }
00094 j0 = e1*znm1+((double)n-1.0e0)*j0;
00095 j1 = e1*zn+(double)n*j1;
00096 znm1 = z2*znm1;
00097 zn = z2*zn;
00098 w = w0*w;
00099 t0 = d[n-1]*w*j0;
00100 w = w0*w;
00101 t1 = d[np1-1]*w*j1;
00102 sum += (t0+t1);
00103 if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
00104 }
00105 S80:
00106 u = exp(-bcorr(a,b));
00107 basym = e0*t*u*sum;
00108 return basym;
00109 } /* END */
|
