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