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 */ |