Doxygen Source Code Documentation
ibeta.c File Reference
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
Go to the source code of this file.
Defines | |
#define | NJ 9 |
#define | NPT ((1<<NJ)*10) |
#define | AITKEN(x, y, z) ((x) - ((y)-(x))*((y)-(x)) / (((z)-(y))-((y)-(x))) ) |
#define | DMP(tbl, jt) |
Functions | |
void | ibeta (double xcut, double a, double b, double *f1, double *flogx, double *flog1x) |
int | main (int argc, char *argv[]) |
Define Documentation
|
|
|
Value: do{ fprintf(stderr," table " #tbl ":") ; \ for( jj=0 ; jj <= jt ; jj++ ) fprintf(stderr," %12.5g",tbl[jj]) ; \ fprintf(stderr,"\n") ; } while(0) |
|
Definition at line 6 of file ibeta.c. Referenced by ibeta(). |
|
Definition at line 9 of file ibeta.c. Referenced by ibeta(). |
Function Documentation
|
Definition at line 11 of file ibeta.c. Referenced by main(), and spmpar().
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 /* fprintf(stderr,"ii=%d dx=%g a1=%g b1=%g x=%g\n",ii,dx,a1,b1,x) ; */ 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 /* fprintf(stderr,"ii=%d x=%g xab=%g xabx=%g xab1x=%g a1=%g b1=%g dx=%g\n", 00034 ii,x,xab[ii],xabx[ii],xab1x[ii],a1,b1,dx) ; */ 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 } |
|
\** File : SUMA.c
Input paramters :
Definition at line 78 of file ibeta.c. References a, argc, ibeta(), and strtod().
|