Revision 692 branches/merge-trunk-davef/src/df1b2-separable/df1b2f20.cpp

df1b2f20.cpp (revision 692)
9 9
 * Description not yet available.
10 10
 */
11 11
#include <df1b2fun.h>
12
#include "df3fun.h"
12 13

  
14
df3_one_variable lgam(const df3_one_variable& _x);
15

  
13 16
/**
14
 * Description not yet available.
15
 * \param
16
 */
17
df1b2variable gammlnguts(const df1b2variable _zz)
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)
18 25
{
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.9843695780195716e-6, 
37
    1.5056327351493116e-7};
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++)
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)
45 33
  {
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);
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;
65 44

  
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++)
45
  }
46
  else if (value(v1)==2.0)
70 47
  {
71
    *zd++ =dfx * *xd++;
48
    // value of lgam(2.0) is 0 and
49
    // 1st derivative is 1-phi
50
    // 2nd deriv is PI*PI/6 - 1
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff