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