Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
ibeta.c
Go to the documentation of this file.00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <math.h>
00004
00005 #undef NJ
00006 #define NJ 9
00007
00008 #undef NPT
00009 #define NPT ((1<<NJ)*10)
00010
00011 void ibeta( double xcut , double a , double b ,
00012 double *f1 , double *flogx , double *flog1x )
00013 {
00014 int nn , ii , jj , istep ;
00015 double dx , s1,sx,s1x , x, a1=a-1.0,b1=b-1.0 ;
00016 double as1[NJ+1] , asx[NJ+1] , as1x[NJ+1] ;
00017 double bs1[NJ-1] , bsx[NJ-1] , bs1x[NJ-1] ;
00018 double cs1[NJ-3] , csx[NJ-3] , cs1x[NJ-3] ;
00019 double xab[NPT] , xabx[NPT] , xab1x[NPT] ;
00020
00021 if( xcut <= 0.0 || xcut >= 1.0 || a <= 0.0 || b <= 0.0 ) return ;
00022
00023 if( f1 == NULL && flogx == NULL && flog1x == NULL ) return ;
00024
00025 dx = xcut / NPT ;
00026 for( ii=1 ; ii <= NPT ; ii++ ){
00027 x = ii*dx ;
00028
00029 xab[ii] = pow(x,a1) * pow(1.0-x,b1) ;
00030 xabx[ii] = xab[ii] * log(x) ;
00031 xab1x[ii] = xab[ii] * log(1.0-x) ;
00032
00033
00034
00035 }
00036
00037 for( nn=NPT,istep=1,jj=NJ ; jj >= 0 ; jj--,istep*=2,nn/=2 ){
00038 s1 = sx = s1x = 0.0 ;
00039 for( ii=istep ; ii <= NPT ; ii+=istep ){
00040 s1 += xab[ii] ; sx += xabx[ii] ; s1x += xab1x[ii] ;
00041 }
00042 dx = xcut / nn ;
00043 as1[jj] = s1*dx ; asx[jj] = sx*dx ; as1x[jj] = s1x*dx ;
00044 }
00045
00046 #undef AITKEN
00047 #define AITKEN(x,y,z) ((x) - ((y)-(x))*((y)-(x)) / (((z)-(y))-((y)-(x))) )
00048
00049 for( jj=0 ; jj <= NJ-2 ; jj++ ){
00050 bs1[jj] = AITKEN( as1[jj] , as1[jj+1] , as1[jj+2] ) ;
00051 bsx[jj] = AITKEN( asx[jj] , asx[jj+1] , asx[jj+2] ) ;
00052 bs1x[jj] = AITKEN( as1x[jj] , as1x[jj+1] , as1x[jj+2] ) ;
00053 }
00054
00055 for( jj=0 ; jj <= NJ-4 ; jj++ ){
00056 cs1[jj] = AITKEN( bs1[jj] , bs1[jj+1] , bs1[jj+2] ) ;
00057 csx[jj] = AITKEN( bsx[jj] , bsx[jj+1] , bsx[jj+2] ) ;
00058 cs1x[jj] = AITKEN( bs1x[jj] , bs1x[jj+1] , bs1x[jj+2] ) ;
00059 }
00060
00061 #define DMP(tbl,jt) \
00062 do{ fprintf(stderr," table " #tbl ":") ; \
00063 for( jj=0 ; jj <= jt ; jj++ ) fprintf(stderr," %12.5g",tbl[jj]) ; \
00064 fprintf(stderr,"\n") ; } while(0)
00065
00066 fprintf(stderr,"ibeta(xc=%g,a=%g,b=%g)\n",xcut,a,b) ;
00067 DMP(as1,NJ) ; DMP(bs1,NJ-2) ; DMP(cs1,NJ-4) ;
00068 DMP(asx,NJ) ; DMP(bsx,NJ-2) ; DMP(csx,NJ-4) ;
00069 DMP(as1x,NJ) ; DMP(bs1x,NJ-2) ; DMP(cs1x,NJ-4) ;
00070
00071 if( f1 != NULL ) *f1 = cs1[NJ-4] ;
00072 if( flogx != NULL ) *flogx = csx[NJ-4] ;
00073 if( flog1x != NULL ) *flog1x = cs1x[NJ-4] ;
00074
00075 return ;
00076 }
00077
00078 int main( int argc , char * argv[] )
00079 {
00080 double xcut , a , b , f1,f2,f3 ;
00081
00082 if( argc < 4 ) exit(0) ;
00083 xcut = strtod(argv[1],NULL) ;
00084 a = strtod(argv[2],NULL) ;
00085 b = strtod(argv[3],NULL) ;
00086
00087 ibeta( xcut,a,b , &f1,&f2,&f3 ) ; exit(0) ;
00088 }