-
Notifications
You must be signed in to change notification settings - Fork 0
/
hermite.C
111 lines (95 loc) · 2.07 KB
/
hermite.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
/*
* hermite.C
*
* Green's function for the quantum harmonic oscillator
* K(x,y)= sum_n H_n(x) H_n(y) / w_n where w_n = n+1/2
*
* Linas Vepstas March 2006
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "brat.h"
/* compute hermite polynomial (physicists norm) up to order
* nmax, using simple recursive algo, and place results in array "vals"
*/
void hermite_array (double x, int nmax, double *H_n)
{
H_n[0] = 0.0;
H_n[1] = 2.0*x;
int n;
for (n=1; n<nmax-1; n++)
{
H_n[n+1] = 2.0*(x*H_n[n] - n*H_n[n-1]);
}
}
double **hermite_grid;
int nsteps;
int nmax;
double xmax;
double xdelta;
void init_hermite_grid (double exmax, double exdelta, int enmax)
{
nmax = enmax;
xmax = exmax;
xdelta = exdelta;
nsteps = (int) (xmax/xdelta) +1;
hermite_grid = (double **) malloc (nsteps * sizeof (double *));
double *arr = (double *) malloc (nsteps *nmax*sizeof (double));
double x = 0.0;
int i;
for (i=0; i<nsteps; i++)
{
hermite_grid[i] = arr;
hermite_array (x, nmax, arr);
arr += nmax;
x += xdelta;
}
}
double hermite (double x, int n)
{
double sgn = 1.0;
if (x < 0.0)
{
x = -x;
if (n%2) sgn = -1.0;
}
int iarr = (int) ((x+0.5*xdelta)/xdelta);
return sgn*hermite_grid[iarr][n];
}
static double hermite_green (double x, double y)
{
double acc = 0.0;
int n;
double en = 0.5;
double tn = 1.0;
double nfac = 1.0;
for (n=0; n<nmax; n++)
{
double term = hermite (x,n) * hermite (y,n);
// term /= tn * nfac * en;
term /= tn * nfac;
acc += term;
nfac *= n+1;
en = n+0.5;
tn *= 2.0;
// printf ("duude n=%d term=%g\n", n, term);
}
acc *= exp (-0.5*(x*x+y*y));
acc /= sqrt (M_PI);
// printf ("x=%g y=%g v=%g\n", x,y, acc);
return acc;
}
static double hermite_green_wrap (double x, double y, int itermax, double param)
{
static int is_init = 0;
if (!is_init)
{
is_init = 1;
init_hermite_grid (4.0, 0.005 ,itermax);
printf ("done with hermite initialization, max=%d\n", itermax);
}
return hermite_green (x,y);
}
DECL_MAKE_HEIGHT(hermite_green_wrap);
/* --------------------------- END OF LIFE ------------------------- */