Revision 601 trunk/src/df1b2separable/df1b2f20.cpp
df1b2f20.cpp (revision 601)  

9  9 
* Description not yet available. 
10  10 
*/ 
11  11 
#include <df1b2fun.h> 
12 
#include "df3fun.h" 

13  12  
14 
df3_one_variable lgam(const df3_one_variable& _x); 

15  
16  13 
/** 
17 
* Log Gamma Function 

18 
* 

19 
* Used to find the Natural log of the gamma function. 

20 
* 

21 
* \param _x the argument 

22 
* 

23 
*/ 

24 
df1b2variable lgam(const df1b2variable& _v1) 

14 
* Description not yet available. 

15 
* \param 

16 
*/ 

17 
df1b2variable gammlnguts(const df1b2variable _zz) 

25  18 
{ 
26 
df1b2variable tmp; 

27 
tmp = 0.0; 

28 
const double phi = 0.5772156649015328606065121; 

29 
const double zeta = 1.2020569031595942853997382; 

30 
init_df3_one_variable v1(_v1); 

31  
32 
if (value(v1)==1.0) 

19 
ADUNCONST(df1b2variable,zz) 

20 
df1b2variable u; 

21 
double z = value(_zz); 

22 
//double zdot=1.0; 

23 
//double z2dot=0.0; 

24 
//double z3dot=0.0; 

25 
//const double lpi =1.1447298858494001741434272; 

26 
//const double pi =3.1415926535897932384626432; 

27 
const double lpp =0.9189385332046727417803297; 

28 
int n=7; 

29 
const double c[9]={0.99999999999980993, 

30 
676.5203681218851, 

31 
1259.1392167224028, 

32 
771.32342877765313, 

33 
176.61502916214059, 

34 
12.507343278686905, 

35 
0.13857109526572012, 

36 
9.9843695780195716e6, 

37 
1.5056327351493116e7}; 

38 
z=1.0; 

39 
double x=c[0]; 

40 
double xdot=0.0; 

41 
double x2dot=0.0; 

42 
double x3dot=0.0; 

43 
int i; 

44 
for (i=1;i<=n+1;i++) 

33  45 
{ 
34 
// value of lgam(1.0) is 0 

35 
// 1st derivative is phi 

36 
// 2nd deriv is PI*PI/6 

37 
// 3rd deriv is 2*zeta 

38 
df3_one_variable v; 

39 
v = 0.0; 

40 
*v.get_udot() = phi; 

41 
*v.get_udot2() = PI*PI/6; 

42 
*v.get_udot3() = 2*zeta; 

43 
tmp=v; 

46 
double zinv=1.0/(z+i); 

47 
x+=c[i]*zinv; 

48 
//xdot=c[i]/square(z+i)*zdot; since zdot=1.0 

49 
xdot=c[i]*square(zinv); 

50 
x2dot+=2.0*c[i]*cube(zinv); 

51 
x3dot=6.0*c[i]*fourth(zinv); 

52 
} 

53 
double t=z+n+0.5; 

54 
//double tdot=zdot; 

55 
//return lpp + (z+0.5)*log(t) t + log(x); 

56 
double ans= lpp + (z+0.5)*log(t) t + log(x); 

57 
//double dfx=zdot*log(t) + (z+0.5)/t*tdot tdot +xdot/x; 

58 
// since tdot=1.0 

59 
// since zdot=1.0 

60 
double dfx=log(t) + (z+0.5)/t 1.0 +xdot/x; 

61 
double d2f=2.0/t (z+0.5)/square(t)+x2dot/x 

62 
square(xdot)/square(x); 

63 
double d3f=3.0/square(t) + 2.0*(z+0.5)/cube(t)+ x3dot/x 

64 
3.0*x2dot*xdot/square(x)+2.0*cube(xdot)/cube(x); 

44  65  
45 
} 

46 
else if (value(v1)==2.0) 

66 
double * xd=zz.get_u_dot(); 

67 
double * zd=u.get_u_dot(); 

68 
*u.get_u()=ans; 

69 
for (i=0;i<df1b2variable::nvar;i++) 

47  70 
{ 
48 
// value of lgam(2.0) is 0 and 

49 
// 1st derivative is 1phi 

50 
// 2nd deriv is PI*PI/6  1 

51 
// 3rd deriv is 2*(1zeta) 
Also available in: Unified diff