1 /**********************************************************************
2   PhiF.c:
3 
4     PhiF.c is a subroutine to calculate the value of a radial function
5     at R.
6 
7   Log of PhiF.c:
8 
9      7/Apr/2004  Released by T.Ozaki
10 
11 ***********************************************************************/
12 
PhiF(double R,double * phi0,double * MRV,int Grid_Num)13 double PhiF(double R, double *phi0, double *MRV, int Grid_Num)
14 {
15   int mp_min,mp_max,m,po;
16   double h1,h2,h3,f1,f2,f3,f4;
17   double g1,g2,x1,x2,y1,y2,f,df;
18   double a,b,rm,y12,y22,result;
19 
20   mp_min = 0;
21   mp_max = Grid_Num - 1;
22   po = 0;
23 
24   if (MRV[Grid_Num-1]<R){
25     result = 0.0;
26     po = 1;
27   }
28   else if (R<MRV[0]){
29 
30     po = 1;
31     m = 4;
32     rm = MRV[m];
33 
34     h1 = MRV[m-1] - MRV[m-2];
35     h2 = MRV[m]   - MRV[m-1];
36     h3 = MRV[m+1] - MRV[m];
37 
38     f1 = phi0[m-2];
39     f2 = phi0[m-1];
40     f3 = phi0[m];
41     f4 = phi0[m+1];
42 
43     g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2);
44     g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3);
45 
46     x1 = rm - MRV[m-1];
47     x2 = rm - MRV[m];
48     y1 = x1/h2;
49     y2 = x2/h2;
50     y12 = y1*y1;
51     y22 = y2*y2;
52 
53     f =  y22*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
54        + y12*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
55 
56     df = 2.0*y2/h2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
57        + y22*(2.0*f2 + h2*g1)/h2
58        + 2.0*y1/h2*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1)
59        - y12*(2.0*f3 - h2*g2)/h2;
60 
61     a = 0.5*df/rm;
62     b = f - a*rm*rm;
63     result = a*R*R + b;
64   }
65 
66   else{
67     do{
68       m = (mp_min + mp_max)/2;
69       if (MRV[m]<R)
70         mp_min = m;
71       else
72         mp_max = m;
73     }
74     while((mp_max-mp_min)!=1);
75     m = mp_max;
76 
77     if (m<2)
78       m = 2;
79     else if (Grid_Num<=m)
80       m = Grid_Num - 2;
81   }
82 
83   /****************************************************
84                  Spline like interpolation
85   ****************************************************/
86 
87   if (po==0){
88 
89     h1 = MRV[m-1] - MRV[m-2];
90     h2 = MRV[m]   - MRV[m-1];
91     h3 = MRV[m+1] - MRV[m];
92 
93     f1 = phi0[m-2];
94     f2 = phi0[m-1];
95     f3 = phi0[m];
96     f4 = phi0[m+1];
97 
98     /****************************************************
99                    Treatment of edge points
100     ****************************************************/
101 
102     if (m==1){
103       h1 = -(h2+h3);
104       f1 = f4;
105     }
106     if (m==(Grid_Num-1)){
107       h3 = -(h1+h2);
108       f4 = f1;
109     }
110 
111     /****************************************************
112                  Calculate the value at R
113     ****************************************************/
114 
115     g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2);
116     g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3);
117 
118     x1 = R - MRV[m-1];
119     x2 = R - MRV[m];
120     y1 = x1/h2;
121     y2 = x2/h2;
122 
123     f =  y2*y2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)
124        + y1*y1*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);
125 
126     result = f;
127   }
128 
129   return result;
130 }
131 
132