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