1 /**********************************************************************
2 Spherical_Bessel.c:
3
4 Spherical_Bessel.c is a subroutine to calculate the spherical
5 Bessel functions and its derivative from 0 to lmax
6
7 Log of Spherical_Bessel.c:
8
9 08/Nov/2005 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include "openmx_common.h"
17
18 #define xmin 0.0
19 #define asize_lmax 22
20
21
Spherical_Bessel(double x,int lmax,double * sb,double * dsb)22 void Spherical_Bessel( double x, int lmax, double *sb, double *dsb )
23 {
24 int m,n,nmax,p;
25 double tsb[asize_lmax+10];
26 double invx,vsb0,vsb1,vsb2,vsbi;
27 double j0,j1,j0p,j1p,sf,tmp,si,co,ix,ix2;
28
29 if (x<0.0){
30 printf("minus x is invalid for Spherical_Bessel\n");
31 exit(0);
32 }
33
34 /*
35 printf("x=%15.12f\n",x);
36 */
37
38 /* find an appropriate nmax */
39
40 nmax = lmax + (int)(1.5*x) + 10;
41
42 if (nmax<30) nmax = 30;
43
44 if (asize_lmax<(lmax+1)){
45 printf("asize_lmax should be larger than %d in Spherical_Bessel.c\n",lmax+1);
46 exit(0);
47 }
48
49 /* if x is larger than xmin */
50
51 if ( xmin < x ){
52
53 double *tmp_array;
54
55 invx = 1.0/x;
56
57 /* allocation of array */
58 tmp_array = (double*)malloc(sizeof(double)*(nmax+1));
59
60 for ( n=0; n<nmax; n++){
61 tmp_array[n] = (2.0*(double)n + 1.0)*invx;
62 }
63
64 /* initial values */
65
66 vsb0 = 0.0;
67 vsb1 = 1.0e-14;
68
69 /* downward recurrence from nmax-2 to lmax+2 */
70
71 for ( n=nmax-1; (lmax+2)<n; n-- ){
72
73 vsb2 = tmp_array[n]*vsb1 - vsb0;
74
75 if (0.0<(vsb2-1.0e+250)){
76 tmp = 1.0/vsb2;
77 vsb1 *= tmp;
78 vsb2 = 1.0;
79 }
80
81 vsbi = vsb0;
82 vsb0 = vsb1;
83 vsb1 = vsb2;
84 }
85
86 /* downward recurrence from lmax+1 to 0 */
87
88 n = lmax + 3;
89 tsb[n-1] = vsb1;
90 tsb[n ] = vsb0;
91 tsb[n+1] = vsbi;
92
93 tmp = tsb[n-1];
94 tsb[n-1] /= tmp;
95 tsb[n ] /= tmp;
96
97 for ( n=lmax+2; 0<n; n-- ){
98
99 tsb[n-1] = tmp_array[n]*tsb[n] - tsb[n+1];
100
101 if (1.0e+250<tsb[n-1]){
102 tmp = tsb[n-1];
103 for (m=n-1; m<=lmax+1; m++){
104 tsb[m] /= tmp;
105 }
106 }
107 }
108
109 /* normalization */
110
111 si = sin(x);
112 co = cos(x);
113 ix = 1.0/x;
114 ix2 = ix*ix;
115 j0 = si*ix;
116 j1 = si*ix*ix - co*ix;
117
118 if (fabs(tsb[1])<fabs(tsb[0])) sf = j0/tsb[0];
119 else sf = j1/tsb[1];
120
121 /* tsb to sb */
122
123 for ( n=0; n<=lmax+1; n++ ){
124 sb[n] = tsb[n]*sf;
125 }
126
127 /* derivative of sb */
128
129 dsb[0] = co*ix - si*ix*ix;
130 for ( n=1; n<=lmax; n++ ){
131 dsb[n] = ( (double)n*sb[n-1] - (double)(n+1.0)*sb[n+1] )/(2.0*(double)n + 1.0);
132 }
133
134 /* freeing of array */
135 free(tmp_array);
136
137 }
138
139 /* if x is smaller than xmin */
140
141 else {
142
143 /* sb */
144
145 for ( n=0; n<=lmax; n++ ){
146 sb[n] = 0.0;
147 }
148 sb[0] = 1.0;
149
150 /* derivative of sb */
151
152 dsb[0] = 0.0;
153 for ( n=1; n<=lmax; n++ ){
154 dsb[n] = ( (double)n*sb[n-1] - (double)(n+1.0)*sb[n+1] )/(2.0*(double)n + 1.0);
155 }
156 }
157 }
158
159