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