1 /**
2 * @file MMCollisionInt.cpp
3 */
4
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at https://cantera.org/license.txt for license and copyright information.
7
8 #include "MMCollisionInt.h"
9 #include "cantera/base/utilities.h"
10 #include "cantera/numerics/polyfit.h"
11 #include "cantera/base/stringUtils.h"
12 #include "cantera/base/global.h"
13
14 using namespace std;
15
16 namespace Cantera
17 {
18 const int DeltaDegree = 6;
19
20 double MMCollisionInt::delta[8] = {0.0, 0.25, 0.50, 0.75, 1.0,
21 1.5, 2.0, 2.5
22 };
23
quadInterp(doublereal x0,doublereal * x,doublereal * y)24 doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y)
25 {
26 double dx21 = x[1] - x[0];
27 double dx32 = x[2] - x[1];
28 double dx31 = dx21 + dx32;
29 double dy32 = y[2] - y[1];
30 double dy21 = y[1] - y[0];
31 double a = (dx21*dy32 - dy21*dx32)/(dx21*dx31*dx32);
32 return a*(x0 - x[0])*(x0 - x[1]) + (dy21/dx21)*(x0 - x[1]) + y[1];
33 }
34
35 double MMCollisionInt::tstar22[37] = {
36 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
37 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0,
38 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0,
39 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 75.0, 100.0
40 };
41
42 double MMCollisionInt::omega22_table[37*8] = {
43 4.1005, 4.266, 4.833, 5.742, 6.729, 8.624, 10.34, 11.89,
44 3.2626, 3.305, 3.516, 3.914, 4.433, 5.57, 6.637, 7.618,
45 2.8399, 2.836, 2.936, 3.168, 3.511, 4.329, 5.126, 5.874,
46 2.531, 2.522, 2.586, 2.749, 3.004, 3.64, 4.282, 4.895,
47 2.2837, 2.277, 2.329, 2.46, 2.665, 3.187, 3.727, 4.249,
48 2.0838, 2.081, 2.13, 2.243, 2.417, 2.862, 3.329, 3.786,
49 1.922, 1.924, 1.97, 2.072, 2.225, 2.614, 3.028, 3.435,
50 1.7902, 1.795, 1.84, 1.934, 2.07, 2.417, 2.788, 3.156,
51 1.6823, 1.689, 1.733, 1.82, 1.944, 2.258, 2.596, 2.933,
52 1.5929, 1.601, 1.644, 1.725, 1.838, 2.124, 2.435, 2.746,
53 1.4551, 1.465, 1.504, 1.574, 1.67, 1.913, 2.181, 2.451,
54 1.3551, 1.365, 1.4, 1.461, 1.544, 1.754, 1.989, 2.228,
55 1.28, 1.289, 1.321, 1.374, 1.447, 1.63, 1.838, 2.053,
56 1.2219, 1.231, 1.259, 1.306, 1.37, 1.532, 1.718, 1.912,
57 1.1757, 1.184, 1.209, 1.251, 1.307, 1.451, 1.618, 1.795,
58 1.0933, 1.1, 1.119, 1.15, 1.193, 1.304, 1.435, 1.578,
59 1.0388, 1.044, 1.059, 1.083, 1.117, 1.204, 1.31, 1.428,
60 0.99963, 1.004, 1.016, 1.035, 1.062, 1.133, 1.22, 1.319,
61 0.96988, 0.9732, 0.983, 0.9991, 1.021, 1.079, 1.153, 1.236,
62 0.92676, 0.9291, 0.936, 0.9473, 0.9628, 1.005, 1.058, 1.121,
63 0.89616, 0.8979, 0.903, 0.9114, 0.923, 0.9545, 0.9955, 1.044,
64 0.87272, 0.8741, 0.878, 0.8845, 0.8935, 0.9181, 0.9505, 0.9893,
65 0.85379, 0.8549, 0.858, 0.8632, 0.8703, 0.8901, 0.9164, 0.9482,
66 0.83795, 0.8388, 0.8414, 0.8456, 0.8515, 0.8678, 0.8895, 0.916,
67 0.82435, 0.8251, 0.8273, 0.8308, 0.8356, 0.8493, 0.8676, 0.8901,
68 0.80184, 0.8024, 0.8039, 0.8065, 0.8101, 0.8201, 0.8337, 0.8504,
69 0.78363, 0.784, 0.7852, 0.7872, 0.7899, 0.7976, 0.8081, 0.8212,
70 0.76834, 0.7687, 0.7696, 0.7712, 0.7733, 0.7794, 0.7878, 0.7983,
71 0.75518, 0.7554, 0.7562, 0.7575, 0.7592, 0.7642, 0.7711, 0.7797,
72 0.74364, 0.7438, 0.7445, 0.7455, 0.747, 0.7512, 0.7569, 0.7642,
73 0.71982, 0.72, 0.7204, 0.7211, 0.7221, 0.725, 0.7289, 0.7339,
74 0.70097, 0.7011, 0.7014, 0.7019, 0.7026, 0.7047, 0.7076, 0.7112,
75 0.68545, 0.6855, 0.6858, 0.6861, 0.6867, 0.6883, 0.6905, 0.6932,
76 0.67232, 0.6724, 0.6726, 0.6728, 0.6733, 0.6743, 0.6762, 0.6784,
77 0.65099, 0.651, 0.6512, 0.6513, 0.6516, 0.6524, 0.6534, 0.6546,
78 0.61397, 0.6141, 0.6143, 0.6145, 0.6147, 0.6148, 0.6148, 0.6147,
79 0.5887, 0.5889, 0.5894, 0.59, 0.5903, 0.5901, 0.5895, 0.5885
80 };
81
82 // changed upper limit to 500 from 1.0e10 dgg 5/21/04
83 double MMCollisionInt::tstar[39] = {
84 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
85 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0,
86 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0,
87 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 75.0, 100.0, 500.0
88 };
89
90 double MMCollisionInt::astar_table[39*8] = {
91 1.0065, 1.0840, 1.0840, 1.0840, 1.0840, 1.0840, 1.0840, 1.0840,
92 1.0231, 1.0660, 1.0380, 1.0400, 1.0430, 1.0500, 1.0520, 1.0510,
93 1.0424, 1.0450, 1.0480, 1.0520, 1.0560, 1.0650, 1.0660, 1.0640,
94 1.0719, 1.0670, 1.0600, 1.0550, 1.0580, 1.0680, 1.0710, 1.0710,
95 1.0936, 1.0870, 1.0770, 1.0690, 1.0680, 1.0750, 1.0780, 1.0780,
96 1.1053, 1.0980, 1.0880, 1.0800, 1.0780, 1.0820, 1.0840, 1.0840,
97 1.1104, 1.1040, 1.0960, 1.0890, 1.0860, 1.0890, 1.0900, 1.0900,
98 1.1114, 1.1070, 1.1000, 1.0950, 1.0930, 1.0950, 1.0960, 1.0950,
99 1.1104, 1.1070, 1.1020, 1.0990, 1.0980, 1.1000, 1.1000, 1.0990,
100 1.1086, 1.1060, 1.1020, 1.1010, 1.1010, 1.1050, 1.1050, 1.1040,
101 1.1063, 1.1040, 1.1030, 1.1030, 1.1040, 1.1080, 1.1090, 1.1080,
102 1.1020, 1.1020, 1.1030, 1.1050, 1.1070, 1.1120, 1.1150, 1.1150,
103 1.0985, 1.0990, 1.1010, 1.1040, 1.1080, 1.1150, 1.1190, 1.1200,
104 1.0960, 1.0960, 1.0990, 1.1030, 1.1080, 1.1160, 1.1210, 1.1240,
105 1.0943, 1.0950, 1.0990, 1.1020, 1.1080, 1.1170, 1.1230, 1.1260,
106 1.0934, 1.0940, 1.0970, 1.1020, 1.1070, 1.1160, 1.1230, 1.1280,
107 1.0926, 1.0940, 1.0970, 1.0990, 1.1050, 1.1150, 1.1230, 1.1300,
108 1.0934, 1.0950, 1.0970, 1.0990, 1.1040, 1.1130, 1.1220, 1.1290,
109 1.0948, 1.0960, 1.0980, 1.1000, 1.1030, 1.1120, 1.1190, 1.1270,
110 1.0965, 1.0970, 1.0990, 1.1010, 1.1040, 1.1100, 1.1180, 1.1260,
111 1.0997, 1.1000, 1.1010, 1.1020, 1.1050, 1.1100, 1.1160, 1.1230,
112 1.1025, 1.1030, 1.1040, 1.1050, 1.1060, 1.1100, 1.1150, 1.1210,
113 1.1050, 1.1050, 1.1060, 1.1070, 1.1080, 1.1110, 1.1150, 1.1200,
114 1.1072, 1.1070, 1.1080, 1.1080, 1.1090, 1.1120, 1.1150, 1.1190,
115 1.1091, 1.1090, 1.1090, 1.1100, 1.1110, 1.1130, 1.1150, 1.1190,
116 1.1107, 1.1110, 1.1110, 1.1110, 1.1120, 1.1140, 1.1160, 1.1190,
117 1.1133, 1.1140, 1.1130, 1.1140, 1.1140, 1.1150, 1.1170, 1.1190,
118 1.1154, 1.1150, 1.1160, 1.1160, 1.1160, 1.1170, 1.1180, 1.1200,
119 1.1172, 1.1170, 1.1170, 1.1180, 1.1180, 1.1180, 1.1190, 1.1200,
120 1.1186, 1.1190, 1.1190, 1.1190, 1.1190, 1.1190, 1.1200, 1.1210,
121 1.1199, 1.1200, 1.1200, 1.1200, 1.1200, 1.1210, 1.1210, 1.1220,
122 1.1223, 1.1220, 1.1220, 1.1220, 1.1220, 1.1230, 1.1230, 1.1240,
123 1.1243, 1.1240, 1.1240, 1.1240, 1.1240, 1.1240, 1.1250, 1.1250,
124 1.1259, 1.1260, 1.1260, 1.1260, 1.1260, 1.1260, 1.1260, 1.1260,
125 1.1273, 1.1270, 1.1270, 1.1270, 1.1270, 1.1270, 1.1270, 1.1280,
126 1.1297, 1.1300, 1.1300, 1.1300, 1.1300, 1.1300, 1.1300, 1.1290,
127 1.1339, 1.1340, 1.1340, 1.1350, 1.1350, 1.1340, 1.1340, 1.1320,
128 1.1364, 1.1370, 1.1370, 1.1380, 1.1390, 1.1380, 1.1370, 1.1350,
129 1.14187, 1.14187, 1.14187, 1.14187, 1.14187, 1.14187, 1.14187,
130 1.14187
131 };
132
133 double MMCollisionInt::bstar_table[39*8] = {
134 1.1852, 1.2963, 1.2963, 1.2963, 1.2963, 1.2963,1.2963, 1.2963,
135 1.1960, 1.216, 1.237, 1.269, 1.285, 1.290, 1.297, 1.294,
136 1.2451, 1.257, 1.340, 1.389, 1.366, 1.327, 1.314, 1.278,
137 1.2900, 1.294, 1.272, 1.258, 1.262, 1.282, 1.290, 1.299,
138 1.2986, 1.291, 1.284, 1.278, 1.277, 1.288, 1.294, 1.297,
139 1.2865, 1.281, 1.276, 1.272, 1.277, 1.286, 1.292, 1.298,
140 1.2665, 1.264, 1.261, 1.263, 1.269, 1.284, 1.292, 1.298,
141 1.2455, 1.244, 1.248, 1.255, 1.262, 1.278, 1.289, 1.296,
142 1.2253, 1.225, 1.234, 1.240, 1.252, 1.271, 1.284, 1.295,
143 1.2078, 1.210, 1.216, 1.227, 1.242, 1.264, 1.281, 1.292,
144 1.1919, 1.192, 1.205, 1.216, 1.230, 1.256, 1.273, 1.287,
145 1.1678, 1.172, 1.181, 1.195, 1.209, 1.237, 1.261, 1.277,
146 1.1496, 1.155, 1.161, 1.174, 1.189, 1.221, 1.246, 1.266,
147 1.1366, 1.141, 1.147, 1.159, 1.174, 1.202, 1.231, 1.256,
148 1.1270, 1.130, 1.138, 1.148, 1.162, 1.191, 1.218, 1.242,
149 1.1197, 1.122, 1.129, 1.140, 1.149, 1.178, 1.205, 1.231,
150 1.1080, 1.110, 1.116, 1.122, 1.132, 1.154, 1.180, 1.205,
151 1.1016, 1.103, 1.107, 1.112, 1.120, 1.138, 1.160, 1.183,
152 1.0980, 1.099, 1.102, 1.106, 1.112, 1.127, 1.145, 1.165,
153 1.0958, 1.097, 1.099, 1.102, 1.107, 1.119, 1.135, 1.153,
154 1.0935, 1.094, 1.095, 1.097, 1.100, 1.109, 1.120, 1.134,
155 1.0925, 1.092, 1.094, 1.095, 1.098, 1.104, 1.112, 1.122,
156 1.0922, 1.092, 1.093, 1.094, 1.096, 1.100, 1.106, 1.115,
157 1.0922, 1.092, 1.093, 1.093, 1.095, 1.098, 1.103, 1.110,
158 1.0923, 1.092, 1.093, 1.093, 1.094, 1.097, 1.101, 1.106,
159 1.0923, 1.092, 1.092, 1.093, 1.094, 1.096, 1.099, 1.103,
160 1.0927, 1.093, 1.093, 1.093, 1.094, 1.095, 1.098, 1.101,
161 1.0930, 1.093, 1.093, 1.093, 1.094, 1.094, 1.096, 1.099,
162 1.0933, 1.094, 1.093, 1.094, 1.094, 1.095, 1.096, 1.098,
163 1.0937, 1.093, 1.094, 1.094, 1.094, 1.094, 1.096, 1.097,
164 1.0939, 1.094, 1.094, 1.094, 1.094, 1.095, 1.095, 1.097,
165 1.0943, 1.094, 1.094, 1.094, 1.095, 1.095, 1.096, 1.096,
166 1.0944, 1.095, 1.094, 1.094, 1.094, 1.095, 1.095, 1.096,
167 1.0944, 1.094, 1.095, 1.094, 1.094, 1.095, 1.096, 1.096,
168 1.0943, 1.095, 1.094, 1.094, 1.095, 1.095, 1.095, 1.095,
169 1.0941, 1.094, 1.094, 1.094, 1.094, 1.094, 1.094, 1.096,
170 1.0947, 1.095, 1.094, 1.094, 1.093, 1.093, 1.094, 1.095,
171 1.0957, 1.095, 1.094, 1.093, 1.092, 1.093, 1.093, 1.094,
172 1.10185, 1.10185, 1.10185, 1.10185, 1.10185, 1.10185, 1.10185,
173 1.10185
174 };
175
176 double MMCollisionInt::cstar_table[39*8] = {
177 0.8889, 0.77778, 0.77778,0.77778,0.77778,0.77778,0.77778,0.77778,
178 0.88575, 0.8988, 0.8378, 0.8029, 0.7876, 0.7805, 0.7799, 0.7801,
179 0.87268, 0.8692,0.8647,0.8479,0.8237,0.7975,0.7881,0.7784,
180 0.85182, 0.8525,0.8366,0.8198,0.8054,0.7903,0.7839,0.782,
181 0.83542, 0.8362,0.8306,0.8196,0.8076,0.7918,0.7842,0.7806,
182 0.82629, 0.8278,0.8252,0.8169,0.8074,0.7916,0.7838,0.7802,
183 0.82299, 0.8249,0.823,0.8165,0.8072,0.7922,0.7839,0.7798,
184 0.82357, 0.8257,0.8241,0.8178,0.8084,0.7927,0.7839,0.7794,
185 0.82657, 0.828,0.8264,0.8199,0.8107,0.7939,0.7842,0.7796,
186 0.8311, 0.8234,0.8295,0.8228,0.8136,0.796,0.7854,0.7798,
187 0.8363, 0.8366,0.8342,0.8267,0.8168,0.7986,0.7864,0.7805,
188 0.84762, 0.8474,0.8438,0.8358,0.825,0.8041,0.7904,0.7822,
189 0.85846, 0.8583,0.853,0.8444,0.8336,0.8118,0.7957,0.7854,
190 0.8684, 0.8674,0.8619,0.8531,0.8423,0.8186,0.8011,0.7898,
191 0.87713, 0.8755,0.8709,0.8616,0.8504,0.8265,0.8072,0.7939,
192 0.88479, 0.8831,0.8779,0.8695,0.8578,0.8338,0.8133,0.799,
193 0.89972, 0.8986,0.8936,0.8846,0.8742,0.8504,0.8294,0.8125,
194 0.91028, 0.9089,0.9043,0.8967,0.8869,0.8649,0.8438,0.8253,
195 0.91793, 0.9166,0.9125,0.9058,0.897,0.8768,0.8557,0.8372,
196 0.92371, 0.9226,0.9189,0.9128,0.905,0.8861,0.8664,0.8484,
197 0.93135, 0.9304,0.9274,0.9226,0.9164,0.9006,0.8833,0.8662,
198 0.93607, 0.9353,0.9329,0.9291,0.924,0.9109,0.8958,0.8802,
199 0.93927, 0.9387,0.9366,0.9334,0.9292,0.9162,0.905,0.8911,
200 0.94149, 0.9409,0.9393,0.9366,0.9331,0.9236,0.9122,0.8997,
201 0.94306, 0.9426,0.9412,0.9388,0.9357,0.9276,0.9175,0.9065,
202 0.94419, 0.9437,0.9425,0.9406,0.938,0.9308,0.9219,0.9119,
203 0.94571, 0.9455,0.9445,0.943,0.9409,0.9353,0.9283,0.9201,
204 0.94662, 0.9464,0.9456,0.9444,0.9428,0.9382,0.9325,0.9258,
205 0.94723, 0.9471,0.9464,0.9455,0.9442,0.9405,0.9355,0.9298,
206 0.94764, 0.9474,0.9469,0.9462,0.945,0.9418,0.9378,0.9328,
207 0.9479, 0.9478,0.9474,0.9465,0.9457,0.943,0.9394,0.9352,
208 0.94827, 0.9481,0.948,0.9472,0.9467,0.9447,0.9422,0.9391,
209 0.94842, 0.9484,0.9481,0.9478,0.9472,0.9458,0.9437,0.9415,
210 0.94852, 0.9484,0.9483,0.948,0.9475,0.9465,0.9449,0.943,
211 0.94861, 0.9487,0.9484,0.9481,0.9479,0.9468,0.9455,0.943,
212 0.94872, 0.9486,0.9486,0.9483,0.9482,0.9475,0.9464,0.9452,
213 0.94881, 0.9488,0.9489,0.949,0.9487,0.9482,0.9476,0.9468,
214 0.94863, 0.9487,0.9489,0.9491,0.9493,0.9491,0.9483,0.9476,
215 0.94444, 0.94444,0.94444,0.94444,0.94444,0.94444,0.94444,0.94444
216 };
217
init(doublereal tsmin,doublereal tsmax,int log_level)218 void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
219 {
220 m_loglevel = log_level;
221 debuglog("Collision Integral Polynomial Fits\n", m_loglevel > 0);
222 m_nmin = -1;
223 m_nmax = -1;
224
225 for (int n = 0; n < 37; n++) {
226 if (tsmin > tstar[n+1]) {
227 m_nmin = n;
228 }
229 if (tsmax > tstar[n+1]) {
230 m_nmax = n+1;
231 }
232 }
233 if (m_nmin < 0 || m_nmin >= 36 || m_nmax < 0 || m_nmax > 36) {
234 m_nmin = 0;
235 m_nmax = 36;
236 }
237 if (m_loglevel > 0) {
238 writelogf("T*_min = %g\n", tstar[m_nmin + 1]);
239 writelogf("T*_max = %g\n", tstar[m_nmax + 1]);
240 }
241 m_logTemp.resize(37);
242 doublereal e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0;
243
244 if (m_loglevel > 0) {
245 writelog("Collision integral fits at each tabulated T* vs. delta*.\n"
246 "These polynomial fits are used to interpolate between "
247 "columns (delta*)\n in the Monchick and Mason tables."
248 " They are only used for nonzero delta*.\n");
249 if (log_level < 4) {
250 writelog("polynomial coefficients not printed (log_level < 4)\n");
251 }
252 }
253
254 for (int i = 0; i < 37; i++) {
255 m_logTemp[i] = log(tstar[i+1]);
256 vector_fp c(DeltaDegree+1);
257
258 double rmserr = fitDelta(0, i, DeltaDegree, c.data());
259 if (log_level > 3) {
260 writelogf("\ndelta* fit at T* = %.6g\n", tstar[i+1]);
261 writelog("omega22 = [" + vec2str(c) + "]\n");
262 }
263 m_o22poly.push_back(c);
264 e22 = std::max(e22, rmserr);
265
266 rmserr = fitDelta(1, i, DeltaDegree, c.data());
267 m_apoly.push_back(c);
268 if (log_level > 3) {
269 writelog("A* = [" + vec2str(c) + "]\n");
270 }
271 ea = std::max(ea, rmserr);
272
273 rmserr = fitDelta(2, i, DeltaDegree, c.data());
274 m_bpoly.push_back(c);
275 if (log_level > 3) {
276 writelog("B* = [" + vec2str(c) + "]\n");
277 }
278 eb = std::max(eb, rmserr);
279
280 rmserr = fitDelta(3, i, DeltaDegree, c.data());
281 m_cpoly.push_back(c);
282 if (log_level > 3) {
283 writelog("C* = [" + vec2str(c) + "]\n");
284 }
285 ec = std::max(ec, rmserr);
286
287 if (log_level > 0) {
288 writelogf("max RMS errors in fits vs. delta*:\n"
289 " omega_22 = %12.6g \n"
290 " A* = %12.6g \n"
291 " B* = %12.6g \n"
292 " C* = %12.6g \n", e22, ea, eb, ec);
293 }
294 }
295 }
296
fitDelta(int table,int ntstar,int degree,doublereal * c)297 doublereal MMCollisionInt::fitDelta(int table, int ntstar, int degree, doublereal* c)
298 {
299 vector_fp w(8);
300 doublereal* begin = 0;
301 switch (table) {
302 case 0:
303 begin = omega22_table + 8*ntstar;
304 break;
305 case 1:
306 begin = astar_table + 8*(ntstar + 1);
307 break;
308 case 2:
309 begin = bstar_table + 8*(ntstar + 1);
310 break;
311 case 3:
312 begin = cstar_table + 8*(ntstar + 1);
313 break;
314 default:
315 return 0.0;
316 }
317 w[0] = -1.0;
318 return polyfit(8, degree, delta, begin, w.data(), c);
319 }
320
omega22(double ts,double deltastar)321 doublereal MMCollisionInt::omega22(double ts, double deltastar)
322 {
323 int i;
324 for (i = 0; i < 37; i++) {
325 if (ts < tstar22[i]) {
326 break;
327 }
328 }
329 int i1 = std::max(i - 1, 0);
330 int i2 = i1+3;
331 if (i2 > 36) {
332 i2 = 36;
333 i1 = i2 - 3;
334 }
335 vector_fp values(3);
336 for (i = i1; i < i2; i++) {
337 if (deltastar == 0.0) {
338 values[i-i1] = omega22_table[8*i];
339 } else {
340 values[i-i1] = poly5(deltastar, m_o22poly[i].data());
341 }
342 }
343 return quadInterp(log(ts), &m_logTemp[i1], values.data());
344 }
345
astar(double ts,double deltastar)346 doublereal MMCollisionInt::astar(double ts, double deltastar)
347 {
348 int i;
349 for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
350 break;
351 }
352 int i1 = std::max(i - 1, 0);
353 int i2 = i1+3;
354 if (i2 > 36) {
355 i2 = 36;
356 i1 = i2 - 3;
357 }
358 vector_fp values(3);
359 for (i = i1; i < i2; i++) {
360 if (deltastar == 0.0) {
361 values[i-i1] = astar_table[8*(i + 1)];
362 } else {
363 values[i-i1] = poly5(deltastar, m_apoly[i].data());
364 }
365 }
366 return quadInterp(log(ts), &m_logTemp[i1], values.data());
367 }
368
bstar(double ts,double deltastar)369 doublereal MMCollisionInt::bstar(double ts, double deltastar)
370 {
371 int i;
372 for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
373 break;
374 }
375 int i1 = std::max(i - 1, 0);
376 int i2 = i1+3;
377 if (i2 > 36) {
378 i2 = 36;
379 i1 = i2 - 3;
380 }
381 vector_fp values(3);
382 for (i = i1; i < i2; i++) {
383 if (deltastar == 0.0) {
384 values[i-i1] = bstar_table[8*(i + 1)];
385 } else {
386 values[i-i1] = poly5(deltastar, m_bpoly[i].data());
387 }
388 }
389 return quadInterp(log(ts), &m_logTemp[i1], values.data());
390 }
391
cstar(double ts,double deltastar)392 doublereal MMCollisionInt::cstar(double ts, double deltastar)
393 {
394 int i;
395 for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
396 break;
397 }
398 int i1 = std::max(i - 1,0);
399 int i2 = i1+3;
400 if (i2 > 36) {
401 i2 = 36;
402 i1 = i2 - 3;
403 }
404 vector_fp values(3);
405 for (i = i1; i < i2; i++) {
406 if (deltastar == 0.0) {
407 values[i-i1] = cstar_table[8*(i + 1)];
408 } else {
409 values[i-i1] = poly5(deltastar, m_cpoly[i].data());
410 }
411 }
412 return quadInterp(log(ts), &m_logTemp[i1], values.data());
413 }
414
fit_omega22(int degree,doublereal deltastar,doublereal * o22)415 void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
416 doublereal* o22)
417 {
418 int n = m_nmax - m_nmin + 1;
419 vector_fp values(n);
420 vector_fp w(n);
421 doublereal* logT = &m_logTemp[m_nmin];
422 for (int i = 0; i < n; i++) {
423 if (deltastar == 0.0) {
424 values[i] = omega22_table[8*(i + m_nmin)];
425 } else {
426 values[i] = poly5(deltastar, m_o22poly[i+m_nmin].data());
427 }
428 }
429 w[0]= -1.0;
430 double rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22);
431 if (m_loglevel > 0 && rmserr > 0.01) {
432 warn_user("MMCollisionInt::fit_omega22",
433 "RMS error = {:12.6g} in omega_22 fit "
434 "with delta* = {:12.6g}", rmserr, deltastar);
435 }
436 }
437
fit(int degree,doublereal deltastar,doublereal * a,doublereal * b,doublereal * c)438 void MMCollisionInt::fit(int degree, doublereal deltastar,
439 doublereal* a, doublereal* b, doublereal* c)
440 {
441 int n = m_nmax - m_nmin + 1;
442 vector_fp values(n);
443 vector_fp w(n);
444 doublereal* logT = &m_logTemp[m_nmin];
445 for (int i = 0; i < n; i++) {
446 if (deltastar == 0.0) {
447 values[i] = astar_table[8*(i + m_nmin + 1)];
448 } else {
449 values[i] = poly5(deltastar, m_apoly[i+m_nmin].data());
450 }
451 }
452 w[0]= -1.0;
453 double rmserr = polyfit(n, degree, logT, values.data(), w.data(), a);
454
455 for (int i = 0; i < n; i++) {
456 if (deltastar == 0.0) {
457 values[i] = bstar_table[8*(i + m_nmin + 1)];
458 } else {
459 values[i] = poly5(deltastar, m_bpoly[i+m_nmin].data());
460 }
461 }
462 w[0]= -1.0;
463 rmserr = polyfit(n, degree, logT, values.data(), w.data(), b);
464
465 for (int i = 0; i < n; i++) {
466 if (deltastar == 0.0) {
467 values[i] = cstar_table[8*(i + m_nmin + 1)];
468 } else {
469 values[i] = poly5(deltastar, m_cpoly[i+m_nmin].data());
470 }
471 }
472 w[0]= -1.0;
473 rmserr = polyfit(n, degree, logT, values.data(), w.data(), c);
474 if (m_loglevel > 2) {
475 writelogf("\nT* fit at delta* = %.6g\n", deltastar);
476
477 writelog("astar = [" + vec2str(vector_fp(a, a+degree+1))+ "]\n");
478 if (rmserr > 0.01) {
479 warn_user("MMCollisionInt::fit",
480 "RMS error = {:12.6g} for A* fit", rmserr);
481 }
482
483 writelog("bstar = [" + vec2str(vector_fp(b, b+degree+1))+ "]\n");
484 if (rmserr > 0.01) {
485 warn_user("MMCollisionInt::fit",
486 "RMS error = {:12.6g} for B* fit", rmserr);
487 }
488
489 writelog("cstar = [" + vec2str(vector_fp(c, c+degree+1))+ "]\n");
490 if (rmserr > 0.01) {
491 warn_user("MMCollisionInt::fit",
492 "RMS error = {:12.6g} for C* fit", rmserr);
493 }
494 }
495 }
496
497 } // namespace
498