Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_04.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 double basym(double *a,double *b,double *lambda,double *eps)
00003
00004
00005
00006
00007
00008
00009
00010
00011 {
00012 static double e0 = 1.12837916709551e0;
00013 static double e1 = .353553390593274e0;
00014 static int num = 20;
00015
00016
00017
00018
00019
00020
00021
00022
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
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 }