1//---------------------------------------------------------------------
2//---------------------------------------------------------------------
3#include <math.h>
4#include "header.h"
5
6#define dmax1(x,y) ((x)>(y)? (x):(y))
7
8void  set_constants() {
9
10//---------------------------------------------------------------------
11//---------------------------------------------------------------------
12
13
14      ce(1,1)  = 2.0e0;
15      ce(1,2)  = 0.0e0;
16      ce(1,3)  = 0.0e0;
17      ce(1,4)  = 4.0e0;
18      ce(1,5)  = 5.0e0;
19      ce(1,6)  = 3.0e0;
20      ce(1,7)  = 0.5e0;
21      ce(1,8)  = 0.02e0;
22      ce(1,9)  = 0.01e0;
23      ce(1,10) = 0.03e0;
24      ce(1,11) = 0.5e0;
25      ce(1,12) = 0.4e0;
26      ce(1,13) = 0.3e0;
27
28      ce(2,1)  = 1.0e0;
29      ce(2,2)  = 0.0e0;
30      ce(2,3)  = 0.0e0;
31      ce(2,4)  = 0.0e0;
32      ce(2,5)  = 1.0e0;
33      ce(2,6)  = 2.0e0;
34      ce(2,7)  = 3.0e0;
35      ce(2,8)  = 0.01e0;
36      ce(2,9)  = 0.03e0;
37      ce(2,10) = 0.02e0;
38      ce(2,11) = 0.4e0;
39      ce(2,12) = 0.3e0;
40      ce(2,13) = 0.5e0;
41
42      ce(3,1)  = 2.0e0;
43      ce(3,2)  = 2.0e0;
44      ce(3,3)  = 0.0e0;
45      ce(3,4)  = 0.0e0;
46      ce(3,5)  = 0.0e0;
47      ce(3,6)  = 2.0e0;
48      ce(3,7)  = 3.0e0;
49      ce(3,8)  = 0.04e0;
50      ce(3,9)  = 0.03e0;
51      ce(3,10) = 0.05e0;
52      ce(3,11) = 0.3e0;
53      ce(3,12) = 0.5e0;
54      ce(3,13) = 0.4e0;
55
56      ce(4,1)  = 2.0e0;
57      ce(4,2)  = 2.0e0;
58      ce(4,3)  = 0.0e0;
59      ce(4,4)  = 0.0e0;
60      ce(4,5)  = 0.0e0;
61      ce(4,6)  = 2.0e0;
62      ce(4,7)  = 3.0e0;
63      ce(4,8)  = 0.03e0;
64      ce(4,9)  = 0.05e0;
65      ce(4,10) = 0.04e0;
66      ce(4,11) = 0.2e0;
67      ce(4,12) = 0.1e0;
68      ce(4,13) = 0.3e0;
69
70      ce(5,1)  = 5.0e0;
71      ce(5,2)  = 4.0e0;
72      ce(5,3)  = 3.0e0;
73      ce(5,4)  = 2.0e0;
74      ce(5,5)  = 0.1e0;
75      ce(5,6)  = 0.4e0;
76      ce(5,7)  = 0.3e0;
77      ce(5,8)  = 0.05e0;
78      ce(5,9)  = 0.04e0;
79      ce(5,10) = 0.03e0;
80      ce(5,11) = 0.1e0;
81      ce(5,12) = 0.3e0;
82      ce(5,13) = 0.2e0;
83
84      c1 = 1.4e0;
85      c2 = 0.4e0;
86      c3 = 0.1e0;
87      c4 = 1.0e0;
88      c5 = 1.4e0;
89
90      bt = sqrt(0.5e0);
91
92      dnxm1 = 1.0e0 / (double)(grid_points(1)-1);
93      dnym1 = 1.0e0 / (double)(grid_points(2)-1);
94      dnzm1 = 1.0e0 / (double)(grid_points(3)-1);
95
96      c1c2 = c1 * c2;
97      c1c5 = c1 * c5;
98      c3c4 = c3 * c4;
99      c1345 = c1c5 * c3c4;
100
101      conz1 = (1.0e0-c1c5);
102
103      tx1 = 1.0e0 / (dnxm1 * dnxm1);
104      tx2 = 1.0e0 / (2.0e0 * dnxm1);
105      tx3 = 1.0e0 / dnxm1;
106
107      ty1 = 1.0e0 / (dnym1 * dnym1);
108      ty2 = 1.0e0 / (2.0e0 * dnym1);
109      ty3 = 1.0e0 / dnym1;
110
111      tz1 = 1.0e0 / (dnzm1 * dnzm1);
112      tz2 = 1.0e0 / (2.0e0 * dnzm1);
113      tz3 = 1.0e0 / dnzm1;
114
115      dx1 = 0.75e0;
116      dx2 = 0.75e0;
117      dx3 = 0.75e0;
118      dx4 = 0.75e0;
119      dx5 = 0.75e0;
120
121      dy1 = 0.75e0;
122      dy2 = 0.75e0;
123      dy3 = 0.75e0;
124      dy4 = 0.75e0;
125      dy5 = 0.75e0;
126
127      dz1 = 1.0e0;
128      dz2 = 1.0e0;
129      dz3 = 1.0e0;
130      dz4 = 1.0e0;
131      dz5 = 1.0e0;
132
133      dxmax = dmax1(dx3, dx4);
134      dymax = dmax1(dy2, dy4);
135      dzmax = dmax1(dz2, dz3);
136
137      dssp = 0.25e0 * dmax1(dx1, dmax1(dy1, dz1) );
138
139      c4dssp = 4.0e0 * dssp;
140      c5dssp = 5.0e0 * dssp;
141
142      dttx1 = dt*tx1;
143      dttx2 = dt*tx2;
144      dtty1 = dt*ty1;
145      dtty2 = dt*ty2;
146      dttz1 = dt*tz1;
147      dttz2 = dt*tz2;
148
149      c2dttx1 = 2.0e0*dttx1;
150      c2dtty1 = 2.0e0*dtty1;
151      c2dttz1 = 2.0e0*dttz1;
152
153      dtdssp = dt*dssp;
154
155      comz1  = dtdssp;
156      comz4  = 4.0e0*dtdssp;
157      comz5  = 5.0e0*dtdssp;
158      comz6  = 6.0e0*dtdssp;
159
160      c3c4tx3 = c3c4*tx3;
161      c3c4ty3 = c3c4*ty3;
162      c3c4tz3 = c3c4*tz3;
163
164      dx1tx1 = dx1*tx1;
165      dx2tx1 = dx2*tx1;
166      dx3tx1 = dx3*tx1;
167      dx4tx1 = dx4*tx1;
168      dx5tx1 = dx5*tx1;
169
170      dy1ty1 = dy1*ty1;
171      dy2ty1 = dy2*ty1;
172      dy3ty1 = dy3*ty1;
173      dy4ty1 = dy4*ty1;
174      dy5ty1 = dy5*ty1;
175
176      dz1tz1 = dz1*tz1;
177      dz2tz1 = dz2*tz1;
178      dz3tz1 = dz3*tz1;
179      dz4tz1 = dz4*tz1;
180      dz5tz1 = dz5*tz1;
181
182      c2iv  = 2.5e0;
183      con43 = 4.0e0/3.0e0;
184      con16 = 1.0e0/6.0e0;
185
186      xxcon1 = c3c4tx3*con43*tx3;
187      xxcon2 = c3c4tx3*tx3;
188      xxcon3 = c3c4tx3*conz1*tx3;
189      xxcon4 = c3c4tx3*con16*tx3;
190      xxcon5 = c3c4tx3*c1c5*tx3;
191
192      yycon1 = c3c4ty3*con43*ty3;
193      yycon2 = c3c4ty3*ty3;
194      yycon3 = c3c4ty3*conz1*ty3;
195      yycon4 = c3c4ty3*con16*ty3;
196      yycon5 = c3c4ty3*c1c5*ty3;
197
198      zzcon1 = c3c4tz3*con43*tz3;
199      zzcon2 = c3c4tz3*tz3;
200      zzcon3 = c3c4tz3*conz1*tz3;
201      zzcon4 = c3c4tz3*con16*tz3;
202      zzcon5 = c3c4tz3*c1c5*tz3;
203
204      return;
205}
206