1 /**********************************************************************
2 Dr_KumoF.c:
3
4 Dr_KumoF.c is a subroutine to calculate the radial derivative of
5 an interpolated value of yv at x using the "Kumogata" interpolation.
6
7 Log of Dr_KumoF.c:
8
9 26/Feb./2012 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <math.h>
15 #include "openmx_common.h"
16
17 #ifdef MAX
18 #undef MAX
19 #endif
20 #define MAX(a,b) ((a)>(b))? (a):(b)
21
22 #ifdef MIN
23 #undef MIN
24 #endif
25 #define MIN(a,b) ((a)<(b))? (a):(b)
26
Dr_KumoF(int N,double x,double r,double * xv,double * rv,double * yv)27 double Dr_KumoF(int N, double x, double r, double *xv, double *rv, double *yv)
28 {
29
30 if (x<xv[0]){
31
32 int m;
33 double rm,h1,h2,h3,f1,f2,f3,f4,a,b,r;
34 double g1,g2,x1,x2,y1,y2,y12,y22,f,df;
35
36 r = exp(x);
37
38 m = 4;
39 rm = rv[m];
40
41 h1 = rv[m-1] - rv[m-2];
42 h2 = rv[m] - rv[m-1];
43 h3 = rv[m+1] - rv[m];
44
45 f1 = yv[m-2];
46 f2 = yv[m-1];
47 f3 = yv[m];
48 f4 = yv[m+1];
49
50 g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2);
51 g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3);
52
53 x1 = rm - rv[m-1];
54 x2 = rm - rv[m];
55 y1 = x1/h2;
56 y2 = x2/h2;
57 y12 = y1*y1;
58 y22 = y2*y2;
59
60 f = y22*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
61 + y12*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
62
63 df = 2.0*y2/h2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
64 + y22*(2.0*f2 + h2*g1)/h2
65 + 2.0*y1/h2*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1)
66 - y12*(2.0*f3 - h2*g2)/h2;
67
68 a = 0.5*df/rm;
69 b = f - a*rm*rm;
70 return 2.0*a*r;
71 }
72
73 else{
74
75 int i;
76 double t,dt,y;
77 double xmin,xmax,tmp;
78
79 xmin = xv[0];
80 xmax = xv[N-1];
81 x = MIN(x,xmax);
82 x = MAX(x,xmin);
83
84 tmp = ((double)N-1.0)/(xmax-xmin);
85 t = (x-xmin)*tmp;
86 i = floor(t);
87 dt = t - (double)i;
88
89 return 0.5*(( 3.0*(yv[i+3]-yv[i]-3.0*(yv[i+2]-yv[i+1]))*dt
90 +2.0*(-yv[i+3]+4.0*yv[i+2]-5.0*yv[i+1]+2.0*yv[i]))*dt
91 +(yv[i+2]-yv[i]))*tmp/r;
92 }
93
94 }
95