OpenWAM
ACT_Sub_DLL.cpp
1 #include"ACT_Sub_DLL.h"
2 #include "Constantes.h"
3 
6 double min(double a, double b) {
7  double sol = 0.;
8  if(a < b)
9  sol = a;
10  else
11  sol = b;
12  return (sol);
13 }
14 
17 void FUNCTION_FOR_INTERPOLATION(double *interpolated, double *time_interpolated, double *CAD_to_interpolate,
18  double *vector_to_interpolate, int size_interpolated, int size_to_interpolate, double speed) {
19 
20  int auxiliar = 0, counter = 0;
21  auxiliar = 0;
22  for(counter = 0; counter < size_interpolated; counter++) {
23  while((CAD_to_interpolate[auxiliar] * 60. / (360. * speed) < time_interpolated[counter])
24  && (counter < size_interpolated) && (auxiliar < size_to_interpolate - 2)) {
25  auxiliar = auxiliar + 1;
26  }
27  if(auxiliar == 0) {
28  interpolated[counter] = vector_to_interpolate[auxiliar];
29  } else {
30 
31  interpolated[counter] = vector_to_interpolate[auxiliar - 1] + (vector_to_interpolate[auxiliar] -
32  vector_to_interpolate[auxiliar - 1]) * (time_interpolated[counter] - CAD_to_interpolate[auxiliar - 1] * 60. /
33  (360. * speed)) / (CAD_to_interpolate[auxiliar] * 60. / (360. * speed) - CAD_to_interpolate[auxiliar - 1] * 60. /
34  (360. * speed));
35 
36  }
37  }
38 }
39 
42 void CALCULUS_OF_VIRTUAL_VELOCITY(double *inj_velocity, double *virtual_velocity, double *dmf, double *time_vector,
43  double rofuel, double dc, double n_holes, double nozzle_d, double D, int size, double PI, double speed, double *EOI_IM,
44  double inj_num, double *SOI_IM, double Piston_D, double DBowl, double CTM, double *CAD, double Kswirl) {
45  int counter = 0, auxiliar = 0;
46  double *virt, *w;
47  int i = 0;
48 
49  virt = (double*) malloc(size * sizeof(double));
50  w = (double*) malloc(size * sizeof(double));
51 
52  for(counter = 0; counter < size; counter++) {
53  virtual_velocity[counter] = -9999.;
54  w[counter] = (pow((Piston_D / DBowl),
55  2) * 0.75 * (DBowl / 2) * 2 * PI * speed / 60 * pow((1 / cosh(CAD[counter] / 100)), 3) + 1) * CTM * Kswirl;
56  }
57 
58  for(counter = 0; counter < size; counter++) {
59 
60  inj_velocity[counter] = dmf[counter] * 4. / (rofuel * n_holes * dc * PI * pow(nozzle_d, 2));
61 
62  for(i = 0; i < inj_num; i++) {
63  if((time_vector[counter] > (SOI_IM[i] / (speed * 6))) && ((time_vector[counter] - (SOI_IM[i] / (speed * 6))) < 0.001)) {
64  // Correccion de inicio de pulso: durante el primer milisegundo.
65  // Comprobado el 6-10-2010 que en 1 ms. la correccion ya es despreciable (0.68%).
66  virt[counter] = inj_velocity[counter] * (0.5274 * exp(-(time_vector[counter] - (SOI_IM[i] /
67  (speed * 6))) / 0.00023) + 1.);
68  if(virt[counter] > inj_velocity[counter]) {
69  inj_velocity[counter] = virt[counter];
70  }
71  }
72  }
73 
74  for(auxiliar = counter; auxiliar < size; auxiliar++) {
75  virt[auxiliar] = inj_velocity[counter] * exp(-(time_vector[auxiliar] - time_vector[counter]) / (D * exp((
76  time_vector[auxiliar] - time_vector[counter]) / 0.006)));
77  if(virt[auxiliar] > virtual_velocity[auxiliar]) {
78  virtual_velocity[auxiliar] = virt[auxiliar];
79  }
80  }
81 
82  /*
83  if((time_vector[counter]>EOI_IM[aux]/(6*speed))&&(virtual_velocity[counter]>virtual_velocity[counter-1])){
84  virtual_velocity[counter]=virtual_velocity[counter-1];
85  }
86  */
87 
88  } // counter
89 
90  free(virt);
91  free(w);
92 }
93 
96 void CALCULUS_OF_ACCUMULATED_INJ_RATE(double *acu_dmf, double *dmf, double *time_vector, int size) {
97 
98  int counter = 0;
99 
100  for(counter = 0; counter < size; counter++) {
101  if(counter == 0) {
102  acu_dmf[counter] = 0.;
103  } else {
104  acu_dmf[counter] = acu_dmf[counter - 1] + ((time_vector[counter] - time_vector[counter - 1]) * (dmf[counter] -
105  (dmf[counter] - dmf[counter - 1]) / 2.));
106  }
107  }
108 
109  if(acu_dmf[size - 1] > 0) {
110  for(counter = 0; counter < size; counter++) {
111  acu_dmf[counter] = acu_dmf[counter] / acu_dmf[size - 1];
112 
113  }
114  }
115 }
116 
119 void STOICHIOMETRY_CONSTANTS(double HC, double *Kst1, double *Kst2, double *Kst3, double *Kst4, double *Kst5,
120  double *Kst6) {
121 
122  *Kst1 = (1 + HC / 4) * 32 / (12 + HC);
123  *Kst2 = pow(*Kst1, -1);
124  *Kst3 = 11 / (8 * (1 + HC / 4));
125  *Kst4 = 9 * HC / (32 * (1 + HC / 4));
126  *Kst5 = 44 / (12 + HC);
127  *Kst6 = 9 * HC / (12 + HC);
128 
129 }
130 
133 void CALCULUS_OF_NUMBER_ELEMENTS(int *num_i_IM, double *time_vector, int size, double speed, double *SOI_IM,
134  double *EOI_IM, int inj_num) {
135 
136  int time_counter = 0, inj_counter = 0;
137 
138  for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
139  num_i_IM[inj_counter] = 0;
140  for(time_counter = 0; time_counter < size; time_counter++) {
141  if((time_vector[time_counter] >= SOI_IM[inj_counter] * 60. / (360. * speed))
142  && (time_vector[time_counter] <= EOI_IM[inj_counter] * 60. / (360. * speed))) {
143  num_i_IM[inj_counter] = num_i_IM[inj_counter] + 1;
144  }
145  }
146  }
147 }
148 
151 void CALCULUS_OF_POI(double **POI_IM, double **mfuel_i_IM, double **mfuel_ij_IM, double *acu_dmf, double *time_vector,
152  int size, double speed, int *num_i_IM, int num_j, double *SOI_IM, double* EOI_IM, int inj_num,
153  stControlElementComposition **elementcontrol) {
154  int counter = 0, inj_counter = 0, aux = 0;
155  int *auxiliar_numi;
156 
157  auxiliar_numi = (int*) malloc(inj_num * sizeof(int*));
158 
159  for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
160  auxiliar_numi[inj_counter] = -1;
161  }
162 
163  for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
164  for(counter = 0; counter < size; counter++) {
165  if(time_vector[counter] >= SOI_IM[inj_counter] * 60. / (360. * speed)
166  && time_vector[counter] <= EOI_IM[inj_counter] * 60. / (360. * speed)) {
167 
168  auxiliar_numi[inj_counter] = auxiliar_numi[inj_counter] + 1;
169  aux = auxiliar_numi[inj_counter];
170  POI_IM[inj_counter][aux] = time_vector[counter];
171 
172  // element selection for the model control
173  if((acu_dmf[counter - 1] < 0.1) && (acu_dmf[counter] >= 0.1)) {
174  elementcontrol[0][0].inj_number = inj_counter;
175  elementcontrol[0][0].num_i = aux;
176  elementcontrol[0][0].inj_number = inj_counter;
177  elementcontrol[1][0].num_i = aux;
178  elementcontrol[1][0].inj_number = inj_counter;
179  elementcontrol[1][0].num_i = aux;
180  elementcontrol[2][0].num_i = aux;
181  elementcontrol[2][0].inj_number = inj_counter;
182  elementcontrol[2][0].num_i = aux;
183  }
184  if((acu_dmf[counter - 1] < 0.3) && (acu_dmf[counter] >= 0.3)) {
185  elementcontrol[3][0].inj_number = inj_counter;
186  elementcontrol[3][0].num_i = aux;
187  elementcontrol[3][0].inj_number = inj_counter;
188  elementcontrol[4][0].num_i = aux;
189  elementcontrol[4][0].inj_number = inj_counter;
190  elementcontrol[4][0].num_i = aux;
191  elementcontrol[5][0].num_i = aux;
192  elementcontrol[5][0].inj_number = inj_counter;
193  elementcontrol[5][0].num_i = aux;
194  }
195  if((acu_dmf[counter - 1] < 0.5) && (acu_dmf[counter] >= 0.5)) {
196  elementcontrol[6][0].inj_number = inj_counter;
197  elementcontrol[6][0].num_i = aux;
198  elementcontrol[6][0].inj_number = inj_counter;
199  elementcontrol[7][0].num_i = aux;
200  elementcontrol[7][0].inj_number = inj_counter;
201  elementcontrol[7][0].num_i = aux;
202  elementcontrol[8][0].num_i = aux;
203  elementcontrol[8][0].inj_number = inj_counter;
204  elementcontrol[8][0].num_i = aux;
205  }
206  if((acu_dmf[counter - 1] < 0.7) && (acu_dmf[counter] >= 0.7)) {
207  elementcontrol[9][0].inj_number = inj_counter;
208  elementcontrol[9][0].num_i = aux;
209  elementcontrol[9][0].inj_number = inj_counter;
210  elementcontrol[10][0].num_i = aux;
211  elementcontrol[10][0].inj_number = inj_counter;
212  elementcontrol[10][0].num_i = aux;
213  elementcontrol[11][0].num_i = aux;
214  elementcontrol[11][0].inj_number = inj_counter;
215  elementcontrol[11][0].num_i = aux;
216  }
217  if((acu_dmf[counter - 1] < 0.9) && (acu_dmf[counter] >= 0.9)) {
218  elementcontrol[12][0].inj_number = inj_counter;
219  elementcontrol[12][0].num_i = aux;
220  elementcontrol[12][0].inj_number = inj_counter;
221  elementcontrol[13][0].num_i = aux;
222  elementcontrol[13][0].inj_number = inj_counter;
223  elementcontrol[13][0].num_i = aux;
224  elementcontrol[14][0].num_i = aux;
225  elementcontrol[14][0].inj_number = inj_counter;
226  elementcontrol[14][0].num_i = aux;
227  }
228 
229  if(aux == 0) {
230  mfuel_i_IM[inj_counter][aux] = ((acu_dmf[counter] + acu_dmf[counter + 1]) / 2.) - acu_dmf[counter - 1];
231  } else if(aux == num_i_IM[inj_counter] - 1) {
232  mfuel_i_IM[inj_counter][aux] = acu_dmf[counter + 1] - ((acu_dmf[counter - 1] + acu_dmf[counter]) / 2.);
233  } else {
234  mfuel_i_IM[inj_counter][aux] = ((acu_dmf[counter] + acu_dmf[counter + 1]) / 2.) - ((
235  acu_dmf[counter - 1] + acu_dmf[counter]) / 2.);
236  }
237 
238  }
239  }
240  }
241 
242  for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
243  for(counter = 0; counter < num_i_IM[inj_counter]; counter++) {
244  mfuel_ij_IM[inj_counter][counter] = mfuel_i_IM[inj_counter][counter] / num_j;
245  }
246  }
247 
248  free(auxiliar_numi);
249 
250 }
251 
254 void CALCULATE_CYCLE(double *roair, double *CAD, double delta_t, double *V_cyl, double VTDC, int counter, double speed,
255  double *p_cyl, double *HRF, double *acu_mf, double *Mbb, double *acu_Mbb, double AFe, double f, double mfuel,
256  double mEGR, double mairIVC, double *T_cyl, double HP, double *Yair, double *Yfuel, double *Yburned, double *U,
257  double *CV, double *H_cooler, double *H, double T_Evaporation_fuel, double inj_fuel_temp, double PRECISION_ITERATION,
258  double *defor, double *Rmixture, double Atmosphere_press, double *Gamma, double PI, double Runiv, double Piston_D,
259  double S, double Crank_L, double Connecting_Rod_L, double E, double Piston_Axis_D, double Piston_Crown_H, double DBowl,
260  double VBowl, double M_Connecting_Rod, double M_P_R_PA, double MW_air, double MW_fuel, double MW_burned,
261  double C_ESteel, double C_Mech_Defor, double CTM, double WC1A, double WC1B, double C2, double C_MBLBY, double Cbb,
262  double TPIS, double TCYL_HEAD, double TCYL, double *Qcylhead, double *Qcyl, double *Qpis) {
263 
264  double QT = 0.; /* Resultant heat due to fuel combustion (-) */
265  double QC = 0.; /* Instantaneous heat */
266  double average_Volume = 0.;
267  /* Average volume between the actual and the previous step */
268  double average_Temperature = 0.;
269  /* Average temperature between the actual and the previous step */
270  double average_Pressure = 0.;
271  /* Average temperature between the actual and the previous step */
272  double DU = 0.; /* Increment of internal energy */
273  double MCYL = 0.; /* Cylinder mass */
274  double AERR = 0.; /* Error in the energetic balance */
275  int niter = 0; /* Number of iterations */
276  double Wi = 0.; /* Instantaneous work */
277  double uf = 0.; /* internal energy of fuel */
278  double ecg = 0.; /* enthalpy of evaporated fuel */
279  double ent_ref = 0.; /* Reference enthalpy */
280  double HFTiny = 0.; /* Enthalpy due to the temperature of the fuel */
281  double a = 0.; /* auxiliar variable */
282  double delta_CAD = 0.; /* Calculus interval CAD */
283  double UANT = 0.; /* UANT acumulates a variable that is calculated as the internal energy in the
284  previous step with the actual masic ratios because of only is considered the
285  increasement of energy due to the rise of the temperature not the rise of the
286  mass */
287 
288  if(counter == 0) {
289  acu_mf[counter] = acu_mf[counter] * mfuel;
290  }
291  if(counter == 1) {
292  acu_mf[counter] = acu_mf[counter] * mfuel;
293  acu_mf[counter - 1] = acu_mf[counter - 1] * mfuel;
294  }
295  if(counter >= 2) {
296  acu_mf[counter] = acu_mf[counter] * mfuel;
297  acu_mf[counter - 1] = acu_mf[counter - 1] * mfuel;
298  }
299 
300  delta_CAD = delta_t * 360. * speed / 60.;
301  V_cyl[counter] = VOLUME(CAD[counter], VTDC, PI, Piston_D, Crank_L, Connecting_Rod_L, E);
302 
303  if(counter < 2) {
304 
305  DEFORMATIONS(&(V_cyl[counter]), &(defor[counter]), p_cyl[counter], CAD[counter], delta_CAD, speed, PI, Piston_D, S,
306  Connecting_Rod_L, E, Piston_Axis_D, Piston_Crown_H, M_Connecting_Rod,
307  M_P_R_PA, C_ESteel, C_Mech_Defor);
308 
309  MASIC_RATIO(&Yair[counter], &Yfuel[counter], &Yburned[counter], &Rmixture[counter], HRF[counter], acu_mf[counter],
310  acu_Mbb[counter], AFe, f, mfuel, mEGR, mairIVC, Runiv, MW_air, MW_fuel,
311  MW_burned);
312 
313  PROPERTIES(&U[counter], &CV[counter], T_cyl[counter], T_cyl[counter], Yair[counter], Yfuel[counter], Yburned[counter]);
314 
315  MCYL = mairIVC + acu_mf[counter] - acu_Mbb[counter];
316  roair[counter] = MCYL / V_cyl[counter];
317  p_cyl[counter] = MCYL * Rmixture[counter] * T_cyl[counter] / V_cyl[counter];
318 
319  } else {
320 
321  *Gamma = (CV[counter - 1] + Rmixture[counter - 1]) / CV[counter - 1];
322 
323  Mbb[counter] = BLOW_BY(p_cyl[counter - 1], T_cyl[counter - 1], Rmixture[counter - 1], delta_CAD, speed, *Gamma,
324  Atmosphere_press, Piston_D, C_MBLBY, Cbb);
325 
326  acu_Mbb[counter] = acu_Mbb[counter - 1] + Mbb[counter];
327 
328  DEFORMATIONS(&(V_cyl[counter]), &(defor[counter]), p_cyl[counter - 1], CAD[counter], delta_CAD, speed, PI, Piston_D, S,
329  Connecting_Rod_L, E, Piston_Axis_D, Piston_Crown_H, M_Connecting_Rod,
330  M_P_R_PA, C_ESteel, C_Mech_Defor);
331 
332  MASIC_RATIO(&Yair[counter], &Yfuel[counter], &Yburned[counter], &Rmixture[counter], HRF[counter - 1], acu_mf[counter],
333  acu_Mbb[counter], AFe, f, mfuel, mEGR, mairIVC, Runiv, MW_air, MW_fuel,
334  MW_burned);
335 
336  PROPERTIES(&UANT, &CV[counter], T_cyl[counter - 1], T_cyl[counter - 1], Yair[counter], Yfuel[counter],
337  Yburned[counter]);
338 
339  MCYL = mairIVC + acu_mf[counter] - acu_Mbb[counter];
340 
341  roair[counter] = MCYL / V_cyl[counter];
342 
343  average_Volume = 0.5 * (V_cyl[counter] + V_cyl[counter - 1]);
344 
345  H_cooler[counter] = HEAT_COOLER(p_cyl[counter - 1], p_cyl[0], T_cyl[counter - 1], T_cyl[0], average_Volume, V_cyl[0],
346  delta_CAD, speed, VTDC, &H[counter], PI, Piston_D, S, DBowl, VBowl, CTM,
347  WC1A, WC1B, C2, TPIS, TCYL_HEAD, TCYL, CAD[counter], Qcylhead, Qcyl, Qpis, counter);
348 
349  QT = mfuel * HP;
350  QC = (HRF[counter - 1] - HRF[counter - 2]) * QT;
351  DU = -(p_cyl[counter - 1] * (V_cyl[counter] - V_cyl[counter - 1])) + QC - H_cooler[counter];
352 
353  T_cyl[counter] = T_cyl[counter - 1] + DU / MCYL / CV[counter];
354 
355  ecg = -21776.6 * Runiv / MW_fuel + Runiv / MW_fuel * (-4.5826 * T_Evaporation_fuel + 0.12428 / 2. * pow(
356  T_Evaporation_fuel, 2.) - 0.00007233 / 3. * pow(T_Evaporation_fuel, 3.) + 0.000000016269 / 4. * pow(
357  T_Evaporation_fuel, 4.) - 26067.28 / T_Evaporation_fuel);
358 
359  ent_ref = (Runiv / MW_fuel) * ((-4.5826 * 298.) + (0.12428 / 2. * pow(298., 2.)) - (0.00007233 / 3. * pow(298.,
360  3.)) + (0.000000016269 / 4. * pow(298., 4.)) - (26067.28 * pow(298., -1.)));
361  ecg = ecg - ent_ref;
362 
363  HFTiny = -1852564 + 2195 * (inj_fuel_temp);
364  niter = 0;
365  AERR = 0;
366 
367  do {
368 
369  T_cyl[counter] = T_cyl[counter] - AERR / CV[counter] / MCYL;
370  p_cyl[counter] = MCYL * Rmixture[counter] * T_cyl[counter] / V_cyl[counter];
371  PROPERTIES(&U[counter], &CV[counter], T_cyl[counter], T_cyl[counter - 1], Yair[counter], Yfuel[counter],
372  Yburned[counter]);
373 
374  average_Temperature = 0.5 * (T_cyl[counter] + T_cyl[counter - 1]);
375  average_Pressure = 0.5 * (p_cyl[counter] + p_cyl[counter - 1]);
376  H_cooler[counter] = HEAT_COOLER(average_Pressure, p_cyl[0], average_Temperature, T_cyl[0], average_Volume, V_cyl[0],
377  delta_CAD, speed, VTDC, &H[counter], PI, Piston_D, S, DBowl, VBowl,
378  CTM, WC1A, WC1B, C2, TPIS, TCYL_HEAD, TCYL, CAD[counter], Qcylhead, Qcyl, Qpis, counter);
379 
380  Wi = -average_Pressure * (V_cyl[counter] - V_cyl[counter - 1]);
381 
382  PROPERTIES_FUEL(&uf, T_cyl[counter]);
383 
384  AERR = (MCYL * (U[counter] - UANT)) - Wi - QC + H_cooler[counter] - ((HFTiny - uf) *
385  (acu_mf[counter] - acu_mf[counter - 1])) + (Rmixture[counter] * T_cyl[counter] * Mbb[counter]);
386  niter = niter + 1;
387 
388  a = fabs(AERR / U[counter]);
389  if(niter > 40) {
390  a = 0.000000000001;
391  }
392  } while(a > PRECISION_ITERATION);
393  }
394 
395  if(counter == 0) {
396  acu_mf[counter] = acu_mf[counter] / mfuel;
397  }
398  if(counter == 1) {
399  acu_mf[counter] = acu_mf[counter] / mfuel;
400  acu_mf[counter - 1] = acu_mf[counter - 1] / mfuel;
401  }
402  if(counter >= 2) {
403  acu_mf[counter] = acu_mf[counter] / mfuel;
404  acu_mf[counter - 1] = acu_mf[counter - 1] / mfuel;
405  }
406 }
407 
409 void CALCULUS_OF_MEAN_VARIABLES(double *p_cyl, double *T_cyl, double *dp_da_cyl, double *CAD, double *pmax,
410  double *Tmax, double *dp_da_max, double *p_exit, double *T_exit, int size) {
411 
412  int counter = 0, iter = 0, n = 0;
413  double *dp_da_aux;
414 
415  dp_da_aux = (double*) malloc(size * sizeof(double));
416  n = 4;
417 
418  for(counter = 0; counter < size; counter++) {
419 
420  if(p_cyl[counter] > *pmax) {
421  *pmax = p_cyl[counter];
422  }
423  if(T_cyl[counter] > *Tmax) {
424  *Tmax = T_cyl[counter];
425  }
426  if(counter > 0) {
427  dp_da_aux[counter] = (p_cyl[counter] - p_cyl[counter - 1]) / (CAD[counter] - CAD[counter - 1]);
428  }
429 
430  if(counter == size - 1) {
431  *p_exit = p_cyl[counter];
432  *T_exit = T_cyl[counter];
433  }
434 
435  }
436 
437  for(counter = 0; counter < size; counter++)
438  dp_da_cyl[counter] = dp_da_aux[counter];
439 
440  /* dp_da_cyl[0]=dp_da_aux[0];
441  dp_da_cyl[size-1]=dp_da_aux[size-1];
442 
443  for(iter=0;iter<n;iter++){
444  for(counter=1;counter<size-1;counter++){
445  dp_da_cyl[counter]=(dp_da_aux[counter-1]+dp_da_aux[counter]+dp_da_aux[counter+1])/3;
446  }
447  } */ // Quito el filtrado de la funcion, pues ya no hace falta.
448  for(counter = 0; counter < size - 1; counter++) {
449  if(dp_da_cyl[counter] > *dp_da_max) {
450  *dp_da_max = dp_da_cyl[counter];
451  }
452  }
453 
454  free(dp_da_aux);
455 
456 }
457 
460 double VOLUME(double CAD, double VTDC, double PI, double Piston_D, double Crank_L, double Connecting_Rod_L, double E) {
461 
462  double V_cyl = 0., AREA = 0., AUX = 0., A = 0.;
463 
464  A = CAD * PI / 180.;
465  AREA = PI * Piston_D * Piston_D / 4.;
466  AUX = (Crank_L * sqrt(pow(1. + 1. / (Crank_L / Connecting_Rod_L), 2) - pow(E / Crank_L,
467  2.)) - Crank_L * (cos(A) + sqrt(pow(1. / (Crank_L / Connecting_Rod_L), 2.) - pow(sin(A) - E / Crank_L, 2.))));
468  V_cyl = (VTDC + AREA * AUX);
469 
470  return V_cyl;
471 }
472 
475 void DEFORMATIONS(double *V_cyl, double *DEFOR, double p_cyl, double CAD, double delta_CAD, double speed, double PI,
476  double Piston_D, double S, double Connecting_Rod_L, double E, double Piston_Axis_D, double Piston_Crown_H,
477  double M_Connecting_Rod, double M_P_R_PA, double C_ESteel, double C_Mech_Defor) {
478 
479  double AVp = 0.; /* Increseament of volume due to the pressure */
480  double AVi = 0.; /* Increasement of volume due to the inertia */
481  double Acel = 0.; /* Linear acceleration of the piston */
482  double Lactual = 0.;
483  /* Distance between the piston and the piston head at actual step */
484  double Lanterior = 0.;
485  /* Distance between the piston and the piston head at previous step */
486  double Lposterior = 0.;
487  /* Distance between the piston and the piston head at following step */
488  double auxalfa = 0.; /* Auxiliar variable */
489  double Msist = 0.;
490  /* Msist is the mass of the piston+piston axis+rings+ 0.33*connecting rod */
491 
492  AVp = (PI * pow(Piston_D, 2.) / 4.) * C_Mech_Defor * (p_cyl / C_ESteel) * pow(Piston_D / Piston_Axis_D,
493  2.) * (Piston_Crown_H + Connecting_Rod_L + S / 2.);
494 
495  auxalfa = (CAD - delta_CAD) * PI / 180.;
496  Lanterior = (sqrt(pow((S / 2.) + Connecting_Rod_L, 2.) - pow(E,
497  2.))) - ((S / 2.) * cos(auxalfa) + sqrt(pow(Connecting_Rod_L, 2.) - pow((S / 2.) * sin(auxalfa) - E, 2.)));
498 
499  auxalfa = (CAD + delta_CAD) * PI / 180.;
500  Lposterior = (sqrt(pow((S / 2.) + Connecting_Rod_L, 2.) - pow(E,
501  2.))) - ((S / 2.) * cos(auxalfa) + sqrt(pow(Connecting_Rod_L, 2.) - pow((S / 2.) * sin(auxalfa) - E, 2.)));
502 
503  auxalfa = (CAD) * PI / 180.;
504  Lactual = (sqrt(pow((S / 2.) + Connecting_Rod_L, 2.) - pow(E,
505  2.))) - ((S / 2.) * cos(auxalfa) + sqrt(pow(Connecting_Rod_L, 2.) - pow((S / 2.) * sin(auxalfa) - E, 2.)));
506 
507  Acel = (Lanterior + Lposterior - Lactual * 2.) / pow((delta_CAD / (6. * speed)), 2.);
508 
509  Msist = 0.33 * M_Connecting_Rod + M_P_R_PA;
510 
511  AVi = Msist * Acel * C_Mech_Defor * 1. / C_ESteel * pow(Piston_D / Piston_Axis_D,
512  2.) * (Piston_Crown_H + Connecting_Rod_L + S / 2.);
513 
514  *DEFOR = (AVp - AVi) * 1000000;
515  *V_cyl = *V_cyl + (AVp - AVi);
516 
517 }
518 
521 void MASIC_RATIO(double *Yair, double *Yfuel, double *Yburned, double *Rmixture, double HRF, double acu_mf,
522  double acu_Mbb, double AFe, double f, double mfuel, double mEGR, double mairIVC, double Runiv, double MW_air,
523  double MW_fuel, double MW_burned) {
524 
525  double Yq2 = 0., Yq0 = 0.;
526 
527  Yq2 = (1. + (AFe)) / (1. + (f) - ((acu_Mbb / 2.) / mfuel) * (1. / (1. + (mEGR) / (mairIVC))));
528  Yq0 = Yq2 * (mEGR) / (mairIVC);
529 
530  *Yburned = (((mfuel + mfuel * AFe) * HRF) + ((mEGR) * Yq2) - ((acu_Mbb / 2.) * Yq0)) /
531  (acu_mf + mairIVC - acu_Mbb / 2.);
532 
533  if(acu_mf >= mfuel * HRF) {
534  *Yfuel = (acu_mf - mfuel * HRF) / (acu_mf + mairIVC - acu_Mbb);
535  } else {
536  *Yfuel = 0;
537  }
538 
539  *Yair = 1. - *Yfuel - *Yburned;
540  *Rmixture = *Yair * Runiv / MW_air + *Yfuel * Runiv / MW_fuel + *Yburned * Runiv / MW_burned;
541 
542 }
543 
546 void PROPERTIES(double *u, double *CV, double T_cyl, double T_cyl_pre, double Yair, double Yfuel, double Yburned) {
547 
548  double cva = 0., cvf = 0., cvq = 0., average_Temperature = 0., ua = 0., uf = 0., uq = 0.;
549 
550  average_Temperature = 0.5 * (T_cyl + T_cyl_pre);
551 
552  cva = (-10.4199 * pow(average_Temperature, 0.5)) + 2522.88 - (67227.1 * pow(average_Temperature,
553  -0.5)) + (917124.4 * pow(average_Temperature, -1.)) - (4174853.6 * pow(average_Temperature, -1.5));
554  cvf = -256.4 + (6.95372 * average_Temperature) - (0.00404715 * pow(average_Temperature,
555  2.)) + (0.000000910259 * pow(average_Temperature, 3.)) + (1458487. * pow(average_Temperature, -2.));
556  cvq = 641.154 + (0.43045 * average_Temperature) - (0.0001125 * pow(average_Temperature,
557  2.)) + (0.000000008979 * pow(average_Temperature, 3.));
558 
559  *CV = cva * Yair + cvf * Yfuel + cvq * Yburned;
560 
561  ua = -4193697.9 - (6.9466 * pow(T_cyl, 1.5)) + (2522.88 * T_cyl) - (134454.16 * pow(T_cyl,
562  0.5)) + (917124.39 * log(T_cyl)) + (8349707.14 * pow(T_cyl, -0.5));
563  uf = -1445686.1 - (256.4 * T_cyl) + (3.47686 * pow(T_cyl, 2.)) - (0.00134905 * pow(T_cyl,
564  3.)) + (0.000000227565 * pow(T_cyl, 4.)) - (1458487. * pow(T_cyl, -1.));
565  uq = -3251495. + (1028.75 * T_cyl) - (0.15377 * pow(T_cyl, 2.)) + (0.000067895 * pow(T_cyl, 3.));
566 
567  *u = ua * Yair + uf * Yfuel + uq * Yburned;
568 }
569 
572 void PROPERTIES_FUEL(double *uf, double T_cyl) {
573 
574  *uf = -1445686.1 - (256.4 * T_cyl) + (3.47686 * pow(T_cyl, 2.)) - (0.00134905 * pow(T_cyl,
575  3.)) + (0.000000227565 * pow(T_cyl, 4.)) - (1458487. * pow(T_cyl, -1.));
576 }
577 
580 double HEAT_COOLER(double p_cyl, double pressureIVC, double T_cyl, double temperatureIVC, double average_Volume,
581  double volumeIVC, double delta_CAD, double speed, double VTDC, double *H, double PI, double Piston_D, double S,
582  double DBowl, double VBowl, double CTM, double WC1A, double WC1B, double C2, double TPIS, double TCYL_HEAD, double TCYL,
583  double CAD, double *Qcylhead, double *Qcyl, double *Qpis, int counter)
584 
585 {
586 
587  double C1 = 0.;
588  double cm = 0.; /* mean piston speed */
589  double PA = 0.; /* pressure in the cylinder without combustion */
590  double comb = 0.; /* Term of combustion */
591  double Cylinder_capacity = 0.; /* Cylinder capacity(m3) */
592  double Piston_area = 0.; /* Piston area */
593  double Cylinder_head_area = 0.; /* piston head area */
594  double H_cooler = 0.; /* Heat hung over to the cooler */
595 
596  cm = 2. * S * speed / 60.;
597  C1 = CALCULATE_C1(cm, CTM, WC1A, WC1B, Piston_D, DBowl, speed, CAD, PI);
598 
599  PA = pressureIVC * pow(volumeIVC / average_Volume, 1.36);
600 
601  Cylinder_capacity = PI * Piston_D * Piston_D * S / 4.;
602 
603  comb = ((Cylinder_capacity * temperatureIVC) / (pressureIVC * volumeIVC)) * (p_cyl - PA);
604 
605  *H = 0.012 * pow(Piston_D, -0.2) * pow(T_cyl, -0.53) * pow(p_cyl, 0.8) * pow((C1 + C2 * comb), 0.8);
606 
607  CALCULATE_AREAS(&Piston_area, &Cylinder_head_area, PI, Piston_D, DBowl, VBowl);
608 
609  H_cooler = *H * (Piston_area * (T_cyl - TPIS) + Cylinder_head_area * (T_cyl - TCYL_HEAD) + 4. * ((
610  average_Volume - VTDC) / Piston_D * (T_cyl - TCYL))) * (delta_CAD) / 6. / speed;
611 
612  Qcylhead[counter] = *H * Cylinder_head_area * (T_cyl - TCYL_HEAD) * (delta_CAD) / 6. / speed;
613  Qcyl[counter] = *H * 4. * ((average_Volume - VTDC) / Piston_D * (T_cyl - TCYL)) * (delta_CAD) / 6. / speed;
614  Qpis[counter] = *H * Piston_area * (T_cyl - TPIS) * (delta_CAD) / 6. / speed;
615 
616  return H_cooler;
617 }
618 
621 double CALCULATE_C1(double cm, double CTM, double WC1A, double WC1B, double Piston_D, double DBowl, double speed,
622  double CAD, double PI) {
623 
624  double C1 = 0.;
625  double cu = 0.;
626  double KCTM = 0.;
627  double ratio_CTM = 0.;
628  double x_alfa = 0.;
629 
630  KCTM = exp(-0.200679 * pow(CTM, 0.431262));
631 
632  ratio_CTM = pow(DBowl / Piston_D, 2) * (1 / KCTM);
633  x_alfa = ratio_CTM + 1 / (pow(cosh(CAD / 100), 40) + ratio_CTM / (1 - ratio_CTM));
634 
635  cu = x_alfa * (2 * PI * speed / 60) * CTM * pow(Piston_D / DBowl, 2) * KCTM * (DBowl / 2);
636 
637  C1 = WC1A * cm + WC1B * cu;
638 
639  return C1;
640 }
641 
644 void CALCULATE_AREAS(double *Piston_area, double *Cylinder_head_area, double PI, double Piston_D, double DBowl,
645  double VBowl) {
646 
647  *Piston_area = PI * pow(Piston_D, 2.) / 4. + VBowl * 4. / DBowl;
648  *Cylinder_head_area = PI * pow(Piston_D, 2.) / 4.;
649 }
650 
653 double BLOW_BY(double p_cyl, double T_cyl, double Rmixture, double delta_CAD, double speed, double Gamma,
654  double Atmosphere_press, double Piston_D, double C_MBLBY, double Cbb) {
655 
656  double C_Z = 0.;
657  double Pressure_up = 0.;
658  double Pressure_down = 0.;
659  double Pressure_critic = 0.;
660  double BBy = 0.;
661 
662  if(p_cyl > Atmosphere_press) {
663  Pressure_up = p_cyl;
664  Pressure_down = Atmosphere_press;
665  } else {
666  Pressure_up = Atmosphere_press;
667  Pressure_down = p_cyl;
668  }
669 
670  Pressure_critic = Pressure_up * pow(2. / (Gamma + 1.), (Gamma / (Gamma - 1.)));
671 
672  if(Pressure_down < Pressure_critic) {
673  Pressure_down = Pressure_critic;
674  }
675 
676  C_Z = ((2. * Gamma) / (Gamma - 1.)) * (pow(Pressure_down / Pressure_up, 2. / Gamma) - pow(Pressure_down / Pressure_up,
677  ((Gamma + 1.) / Gamma)));
678  BBy = Cbb * C_MBLBY * Piston_D * Pressure_up * pow(C_Z / (Rmixture * T_cyl), 0.5);
679 
680  if(Atmosphere_press > p_cyl) {
681  BBy = -BBy;
682  }
683 
684  BBy = BBy * delta_CAD / 6. / speed;
685 
686  return BBy;
687 }
688 
691 void CALCULUS_OF_IMP_HP(double *complete_p_cyl, double *complete_CAD, double *p_cyl, double *V_cyl,
692  double *complete_V_cyl, double *complete_deform, double *WI_HP, double *IMP_HP, int complete_size,
693  int complete_prev_size, double delta_t, double speed, int size, double IVC, double EVO, double VTDC,
694  double Cylinder_capacity, double PI, double Piston_D, double S, double Crank_L, double Connecting_Rod_L, double E,
695  double Piston_Axis_D, double Piston_Crown_H, double M_Connecting_Rod, double M_P_R_PA, double C_ESteel,
696  double C_Mech_Defor, double inlet_pres, double exhaust_pres) {
697 
698  double delta_CAD = 0.;
699  int counter = 0;
700 
701  *IMP_HP = 0;
702  *WI_HP = 0;
703  delta_CAD = delta_t * 360. * speed / 60.;
704  complete_CAD[0] = IVC - complete_prev_size * delta_t * 360. * speed / 60.;
705 
706  for(counter = 1; counter < complete_size; counter++) {
707  complete_CAD[counter] = complete_CAD[counter - 1] + delta_t * 360. * speed / 60.;
708  }
709 
710  for(counter = 0; counter < complete_size; counter++) {
711  if(counter < complete_prev_size) {
712  complete_p_cyl[counter] = ((inlet_pres - p_cyl[0]) / (-180 - IVC)) * (complete_CAD[counter] - (-180)) + inlet_pres;
713  }
714  if((counter >= complete_prev_size) && (counter < complete_prev_size + size)) {
715  complete_p_cyl[counter] = p_cyl[counter - complete_prev_size];
716  }
717  if(counter >= complete_prev_size + size) {
718  complete_p_cyl[counter] = ((exhaust_pres - p_cyl[size - 1]) / (180 - EVO)) * (complete_CAD[counter] -
719  (180)) + exhaust_pres;
720  }
721  }
722 
723  for(counter = 0; counter < complete_size; counter++) {
724  if(counter < complete_prev_size) {
725  complete_V_cyl[counter] = VOLUME(complete_CAD[counter], VTDC, PI, Piston_D, Crank_L, Connecting_Rod_L, E);
726  DEFORMATIONS(&complete_V_cyl[counter], &complete_deform[counter], complete_p_cyl[counter], complete_CAD[counter],
727  delta_CAD, speed, PI, Piston_D, S, Connecting_Rod_L, E, Piston_Axis_D,
728  Piston_Crown_H, M_Connecting_Rod, M_P_R_PA, C_ESteel, C_Mech_Defor);
729  }
730  if((counter >= complete_prev_size) && (counter < complete_prev_size + size)) {
731  complete_V_cyl[counter] = V_cyl[counter - complete_prev_size];
732  }
733  if(counter >= complete_prev_size + size) {
734  complete_V_cyl[counter] = VOLUME(complete_CAD[counter], VTDC, PI, Piston_D, Crank_L, Connecting_Rod_L, E);
735  DEFORMATIONS(&complete_V_cyl[counter], &complete_deform[counter], complete_p_cyl[counter], complete_CAD[counter],
736  delta_CAD, speed, PI, Piston_D, S, Connecting_Rod_L, E, Piston_Axis_D,
737  Piston_Crown_H, M_Connecting_Rod, M_P_R_PA, C_ESteel, C_Mech_Defor);
738  }
739 
740  if(counter > 0) {
741  *WI_HP = *WI_HP + (complete_p_cyl[counter] + complete_p_cyl[counter - 1]) * 0.5 * (complete_V_cyl[counter] -
742  complete_V_cyl[counter - 1]);
743  }
744  }
745  *IMP_HP = *WI_HP / Cylinder_capacity;
746 
747 }
748 
751 void FUNCTION_NOX(double *YNOeq_value, double *KdYNO_value, double **YNOeq, double **KdYNO, double temperature,
752  double mO2, double mtotal) {
753  int i = 0, j = 0;
754  double di = 0., dj = 0.;
755 
756  YNOeq[0][0] = 0.000000;
757  KdYNO[0][0] = 0.000000;
758  YNOeq[0][1] = 0.000034;
759  KdYNO[0][1] = 0.000000;
760  YNOeq[0][2] = 0.000054;
761  KdYNO[0][2] = 0.000000;
762  YNOeq[0][3] = 0.000074;
763  KdYNO[0][3] = 0.000000;
764  YNOeq[0][4] = 0.000000;
765  KdYNO[0][4] = 0.000000;
766  YNOeq[0][5] = 0.000000;
767  KdYNO[0][5] = 0.000000;
768  YNOeq[0][6] = 0.000000;
769  KdYNO[0][6] = 0.000000;
770  YNOeq[0][7] = 0.000000;
771  KdYNO[0][7] = 0.000000;
772  YNOeq[0][8] = 0.000000;
773  KdYNO[0][8] = 0.000000;
774  YNOeq[0][9] = 0.000000;
775  KdYNO[0][9] = 0.000000;
776  YNOeq[0][10] = 0.000000;
777  KdYNO[0][10] = 0.000000;
778  YNOeq[0][11] = 0.000000;
779  KdYNO[0][11] = 0.000000;
780  YNOeq[1][0] = 0.000000;
781  KdYNO[1][0] = 0.000000;
782  YNOeq[1][1] = 0.000068;
783  KdYNO[1][1] = 0.000000;
784  YNOeq[1][2] = 0.000090;
785  KdYNO[1][2] = 0.000000;
786  YNOeq[1][3] = 0.000113;
787  KdYNO[1][3] = 0.000000;
788  YNOeq[1][4] = 0.000135;
789  KdYNO[1][4] = 0.000000;
790  YNOeq[1][5] = 0.000160;
791  KdYNO[1][5] = 0.000000;
792  YNOeq[1][6] = 0.000186;
793  KdYNO[1][6] = 0.000000;
794  YNOeq[1][7] = 0.000210;
795  KdYNO[1][7] = 0.000000;
796  YNOeq[1][8] = 0.000000;
797  KdYNO[1][8] = 0.000000;
798  YNOeq[1][9] = 0.000000;
799  KdYNO[1][9] = 0.000000;
800  YNOeq[1][10] = 0.000000;
801  KdYNO[1][10] = 0.000000;
802  YNOeq[1][11] = 0.000000;
803  KdYNO[1][11] = 0.000000;
804  YNOeq[2][0] = 0.000000;
805  KdYNO[2][0] = 0.000000;
806  YNOeq[2][1] = 0.000080;
807  KdYNO[2][1] = 0.000000;
808  YNOeq[2][2] = 0.000131;
809  KdYNO[2][2] = 0.000000;
810  YNOeq[2][3] = 0.000168;
811  KdYNO[2][3] = 0.000000;
812  YNOeq[2][4] = 0.000197;
813  KdYNO[2][4] = 0.000000;
814  YNOeq[2][5] = 0.000224;
815  KdYNO[2][5] = 0.000000;
816  YNOeq[2][6] = 0.000251;
817  KdYNO[2][6] = 0.000000;
818  YNOeq[2][7] = 0.000277;
819  KdYNO[2][7] = 0.000000;
820  YNOeq[2][8] = 0.000304;
821  KdYNO[2][8] = 0.000000;
822  YNOeq[2][9] = 0.000000;
823  KdYNO[2][9] = 0.000000;
824  YNOeq[2][10] = 0.000000;
825  KdYNO[2][10] = 0.000000;
826  YNOeq[2][11] = 0.000000;
827  KdYNO[2][11] = 0.000000;
828  YNOeq[3][0] = 0.000000;
829  KdYNO[3][0] = 0.000000;
830  YNOeq[3][1] = 0.000056;
831  KdYNO[3][1] = 0.000000;
832  YNOeq[3][2] = 0.000111;
833  KdYNO[3][2] = 0.000000;
834  YNOeq[3][3] = 0.000165;
835  KdYNO[3][3] = 0.000000;
836  YNOeq[3][4] = 0.000220;
837  KdYNO[3][4] = 0.000000;
838  YNOeq[3][5] = 0.000274;
839  KdYNO[3][5] = 0.000000;
840  YNOeq[3][6] = 0.000329;
841  KdYNO[3][6] = 0.000000;
842  YNOeq[3][7] = 0.000366;
843  KdYNO[3][7] = 0.000000;
844  YNOeq[3][8] = 0.000396;
845  KdYNO[3][8] = 0.000000;
846  YNOeq[3][9] = 0.000000;
847  KdYNO[3][9] = 0.000000;
848  YNOeq[3][10] = 0.000000;
849  KdYNO[3][10] = 0.000000;
850  YNOeq[3][11] = 0.000000;
851  KdYNO[3][11] = 0.000000;
852  YNOeq[4][0] = 0.000001;
853  KdYNO[4][0] = 0.000000;
854  YNOeq[4][1] = 0.000150;
855  KdYNO[4][1] = 0.000000;
856  YNOeq[4][2] = 0.000226;
857  KdYNO[4][2] = 0.000000;
858  YNOeq[4][3] = 0.000280;
859  KdYNO[4][3] = 0.000000;
860  YNOeq[4][4] = 0.000324;
861  KdYNO[4][4] = 0.000000;
862  YNOeq[4][5] = 0.000368;
863  KdYNO[4][5] = 0.000000;
864  YNOeq[4][6] = 0.000411;
865  KdYNO[4][6] = 0.000000;
866  YNOeq[4][7] = 0.000455;
867  KdYNO[4][7] = 0.000000;
868  YNOeq[4][8] = 0.000500;
869  KdYNO[4][8] = 0.000000;
870  YNOeq[4][9] = 0.000551;
871  KdYNO[4][9] = 0.000000;
872  YNOeq[4][10] = 0.000000;
873  KdYNO[4][10] = 0.000000;
874  YNOeq[4][11] = 0.000000;
875  KdYNO[4][11] = 0.000000;
876  YNOeq[5][0] = 0.000001;
877  KdYNO[5][0] = 0.000000;
878  YNOeq[5][1] = 0.000089;
879  KdYNO[5][1] = 0.000000;
880  YNOeq[5][2] = 0.000173;
881  KdYNO[5][2] = 0.000000;
882  YNOeq[5][3] = 0.000258;
883  KdYNO[5][3] = 0.000000;
884  YNOeq[5][4] = 0.000342;
885  KdYNO[5][4] = 0.000000;
886  YNOeq[5][5] = 0.000426;
887  KdYNO[5][5] = 0.000000;
888  YNOeq[5][6] = 0.000510;
889  KdYNO[5][6] = 0.000000;
890  YNOeq[5][7] = 0.000590;
891  KdYNO[5][7] = 0.000000;
892  YNOeq[5][8] = 0.000644;
893  KdYNO[5][8] = 0.000000;
894  YNOeq[5][9] = 0.000696;
895  KdYNO[5][9] = 0.000000;
896  YNOeq[5][10] = 0.000000;
897  KdYNO[5][10] = 0.000000;
898  YNOeq[5][11] = 0.000000;
899  KdYNO[5][11] = 0.000000;
900  YNOeq[6][0] = 0.000002;
901  KdYNO[6][0] = 0.000000;
902  YNOeq[6][1] = 0.000235;
903  KdYNO[6][1] = 0.000000;
904  YNOeq[6][2] = 0.000351;
905  KdYNO[6][2] = 0.000000;
906  YNOeq[6][3] = 0.000438;
907  KdYNO[6][3] = 0.000000;
908  YNOeq[6][4] = 0.000511;
909  KdYNO[6][4] = 0.000000;
910  YNOeq[6][5] = 0.000585;
911  KdYNO[6][5] = 0.000000;
912  YNOeq[6][6] = 0.000658;
913  KdYNO[6][6] = 0.000000;
914  YNOeq[6][7] = 0.000727;
915  KdYNO[6][7] = 0.000000;
916  YNOeq[6][8] = 0.000790;
917  KdYNO[6][8] = 0.000000;
918  YNOeq[6][9] = 0.000858;
919  KdYNO[6][9] = 0.000000;
920  YNOeq[6][10] = 0.000931;
921  KdYNO[6][10] = 0.000000;
922  YNOeq[6][11] = 0.000000;
923  KdYNO[6][11] = 0.000000;
924  YNOeq[7][0] = 0.000005;
925  KdYNO[7][0] = 0.000000;
926  YNOeq[7][1] = 0.000136;
927  KdYNO[7][1] = 0.000000;
928  YNOeq[7][2] = 0.000254;
929  KdYNO[7][2] = 0.000000;
930  YNOeq[7][3] = 0.000372;
931  KdYNO[7][3] = 0.000000;
932  YNOeq[7][4] = 0.000489;
933  KdYNO[7][4] = 0.000000;
934  YNOeq[7][5] = 0.000607;
935  KdYNO[7][5] = 0.000000;
936  YNOeq[7][6] = 0.000725;
937  KdYNO[7][6] = 0.000000;
938  YNOeq[7][7] = 0.000843;
939  KdYNO[7][7] = 0.000001;
940  YNOeq[7][8] = 0.000961;
941  KdYNO[7][8] = 0.000001;
942  YNOeq[7][9] = 0.001066;
943  KdYNO[7][9] = 0.000001;
944  YNOeq[7][10] = 0.001140;
945  KdYNO[7][10] = 0.000001;
946  YNOeq[7][11] = 0.000000;
947  KdYNO[7][11] = 0.000000;
948  YNOeq[8][0] = 0.000010;
949  KdYNO[8][0] = 0.000000;
950  YNOeq[8][1] = 0.000352;
951  KdYNO[8][1] = 0.000000;
952  YNOeq[8][2] = 0.000525;
953  KdYNO[8][2] = 0.000000;
954  YNOeq[8][3] = 0.000653;
955  KdYNO[8][3] = 0.000001;
956  YNOeq[8][4] = 0.000762;
957  KdYNO[8][4] = 0.000001;
958  YNOeq[8][5] = 0.000871;
959  KdYNO[8][5] = 0.000001;
960  YNOeq[8][6] = 0.000981;
961  KdYNO[8][6] = 0.000001;
962  YNOeq[8][7] = 0.001081;
963  KdYNO[8][7] = 0.000002;
964  YNOeq[8][8] = 0.001177;
965  KdYNO[8][8] = 0.000002;
966  YNOeq[8][9] = 0.001270;
967  KdYNO[8][9] = 0.000002;
968  YNOeq[8][10] = 0.001362;
969  KdYNO[8][10] = 0.000002;
970  YNOeq[8][11] = 0.001474;
971  KdYNO[8][11] = 0.000003;
972  YNOeq[9][0] = 0.000017;
973  KdYNO[9][0] = 0.000000;
974  YNOeq[9][1] = 0.000210;
975  KdYNO[9][1] = 0.000001;
976  YNOeq[9][2] = 0.000376;
977  KdYNO[9][2] = 0.000001;
978  YNOeq[9][3] = 0.000543;
979  KdYNO[9][3] = 0.000002;
980  YNOeq[9][4] = 0.000709;
981  KdYNO[9][4] = 0.000003;
982  YNOeq[9][5] = 0.000875;
983  KdYNO[9][5] = 0.000003;
984  YNOeq[9][6] = 0.001041;
985  KdYNO[9][6] = 0.000004;
986  YNOeq[9][7] = 0.001208;
987  KdYNO[9][7] = 0.000005;
988  YNOeq[9][8] = 0.001374;
989  KdYNO[9][8] = 0.000005;
990  YNOeq[9][9] = 0.001523;
991  KdYNO[9][9] = 0.000006;
992  YNOeq[9][10] = 0.001629;
993  KdYNO[9][10] = 0.000007;
994  YNOeq[9][11] = 0.001748;
995  KdYNO[9][11] = 0.000009;
996  YNOeq[10][0] = 0.000029;
997  KdYNO[10][0] = 0.000000;
998  YNOeq[10][1] = 0.000508;
999  KdYNO[10][1] = 0.000003;
1000  YNOeq[10][2] = 0.000755;
1001  KdYNO[10][2] = 0.000005;
1002  YNOeq[10][3] = 0.000936;
1003  KdYNO[10][3] = 0.000006;
1004  YNOeq[10][4] = 0.001085;
1005  KdYNO[10][4] = 0.000008;
1006  YNOeq[10][5] = 0.001234;
1007  KdYNO[10][5] = 0.000010;
1008  YNOeq[10][6] = 0.001382;
1009  KdYNO[10][6] = 0.000012;
1010  YNOeq[10][7] = 0.001530;
1011  KdYNO[10][7] = 0.000014;
1012  YNOeq[10][8] = 0.001676;
1013  KdYNO[10][8] = 0.000016;
1014  YNOeq[10][9] = 0.001806;
1015  KdYNO[10][9] = 0.000018;
1016  YNOeq[10][10] = 0.001926;
1017  KdYNO[10][10] = 0.000020;
1018  YNOeq[10][11] = 0.002061;
1019  KdYNO[10][11] = 0.000025;
1020  YNOeq[11][0] = 0.000045;
1021  KdYNO[11][0] = 0.000001;
1022  YNOeq[11][1] = 0.000307;
1023  KdYNO[11][1] = 0.000005;
1024  YNOeq[11][2] = 0.000535;
1025  KdYNO[11][2] = 0.000010;
1026  YNOeq[11][3] = 0.000763;
1027  KdYNO[11][3] = 0.000015;
1028  YNOeq[11][4] = 0.000992;
1029  KdYNO[11][4] = 0.000020;
1030  YNOeq[11][5] = 0.001220;
1031  KdYNO[11][5] = 0.000026;
1032  YNOeq[11][6] = 0.001448;
1033  KdYNO[11][6] = 0.000031;
1034  YNOeq[11][7] = 0.001676;
1035  KdYNO[11][7] = 0.000036;
1036  YNOeq[11][8] = 0.001905;
1037  KdYNO[11][8] = 0.000041;
1038  YNOeq[11][9] = 0.002108;
1039  KdYNO[11][9] = 0.000046;
1040  YNOeq[11][10] = 0.002269;
1041  KdYNO[11][10] = 0.000056;
1042  YNOeq[11][11] = 0.002411;
1043  KdYNO[11][11] = 0.000067;
1044  YNOeq[12][0] = 0.000067;
1045  KdYNO[12][0] = 0.000002;
1046  YNOeq[12][1] = 0.000709;
1047  KdYNO[12][1] = 0.000024;
1048  YNOeq[12][2] = 0.001051;
1049  KdYNO[12][2] = 0.000036;
1050  YNOeq[12][3] = 0.001300;
1051  KdYNO[12][3] = 0.000047;
1052  YNOeq[12][4] = 0.001498;
1053  KdYNO[12][4] = 0.000059;
1054  YNOeq[12][5] = 0.001696;
1055  KdYNO[12][5] = 0.000070;
1056  YNOeq[12][6] = 0.001894;
1057  KdYNO[12][6] = 0.000082;
1058  YNOeq[12][7] = 0.002101;
1059  KdYNO[12][7] = 0.000098;
1060  YNOeq[12][8] = 0.002312;
1061  KdYNO[12][8] = 0.000116;
1062  YNOeq[12][9] = 0.002488;
1063  KdYNO[12][9] = 0.000131;
1064  YNOeq[12][10] = 0.002642;
1065  KdYNO[12][10] = 0.000143;
1066  YNOeq[12][11] = 0.002805;
1067  KdYNO[12][11] = 0.000166;
1068  YNOeq[13][0] = 0.000093;
1069  KdYNO[13][0] = 0.000007;
1070  YNOeq[13][1] = 0.000441;
1071  KdYNO[13][1] = 0.000035;
1072  YNOeq[13][2] = 0.000745;
1073  KdYNO[13][2] = 0.000068;
1074  YNOeq[13][3] = 0.001049;
1075  KdYNO[13][3] = 0.000100;
1076  YNOeq[13][4] = 0.001354;
1077  KdYNO[13][4] = 0.000133;
1078  YNOeq[13][5] = 0.001658;
1079  KdYNO[13][5] = 0.000165;
1080  YNOeq[13][6] = 0.001962;
1081  KdYNO[13][6] = 0.000198;
1082  YNOeq[13][7] = 0.002267;
1083  KdYNO[13][7] = 0.000230;
1084  YNOeq[13][8] = 0.002571;
1085  KdYNO[13][8] = 0.000263;
1086  YNOeq[13][9] = 0.002840;
1087  KdYNO[13][9] = 0.000297;
1088  YNOeq[13][10] = 0.003048;
1089  KdYNO[13][10] = 0.000349;
1090  YNOeq[13][11] = 0.003239;
1091  KdYNO[13][11] = 0.000415;
1092  YNOeq[14][0] = 0.000132;
1093  KdYNO[14][0] = 0.000021;
1094  YNOeq[14][1] = 0.000958;
1095  KdYNO[14][1] = 0.000156;
1096  YNOeq[14][2] = 0.001419;
1097  KdYNO[14][2] = 0.000238;
1098  YNOeq[14][3] = 0.001753;
1099  KdYNO[14][3] = 0.000305;
1100  YNOeq[14][4] = 0.002011;
1101  KdYNO[14][4] = 0.000368;
1102  YNOeq[14][5] = 0.002269;
1103  KdYNO[14][5] = 0.000430;
1104  YNOeq[14][6] = 0.002527;
1105  KdYNO[14][6] = 0.000492;
1106  YNOeq[14][7] = 0.002807;
1107  KdYNO[14][7] = 0.000597;
1108  YNOeq[14][8] = 0.003099;
1109  KdYNO[14][8] = 0.000724;
1110  YNOeq[14][9] = 0.003337;
1111  KdYNO[14][9] = 0.000815;
1112  YNOeq[14][10] = 0.003527;
1113  KdYNO[14][10] = 0.000865;
1114  YNOeq[14][11] = 0.003721;
1115  KdYNO[14][11] = 0.000958;
1116  YNOeq[15][0] = 0.000175;
1117  KdYNO[15][0] = 0.000057;
1118  YNOeq[15][1] = 0.000626;
1119  KdYNO[15][1] = 0.000212;
1120  YNOeq[15][2] = 0.001158;
1121  KdYNO[15][2] = 0.000372;
1122  YNOeq[15][3] = 0.001555;
1123  KdYNO[15][3] = 0.000550;
1124  YNOeq[15][4] = 0.001934;
1125  KdYNO[15][4] = 0.000729;
1126  YNOeq[15][5] = 0.002314;
1127  KdYNO[15][5] = 0.000908;
1128  YNOeq[15][6] = 0.002693;
1129  KdYNO[15][6] = 0.001088;
1130  YNOeq[15][7] = 0.003072;
1131  KdYNO[15][7] = 0.001267;
1132  YNOeq[15][8] = 0.003452;
1133  KdYNO[15][8] = 0.001447;
1134  YNOeq[15][9] = 0.003733;
1135  KdYNO[15][9] = 0.001636;
1136  YNOeq[15][10] = 0.003997;
1137  KdYNO[15][10] = 0.001895;
1138  YNOeq[15][11] = 0.004246;
1139  KdYNO[15][11] = 0.002222;
1140  YNOeq[16][0] = 0.000207;
1141  KdYNO[16][0] = 0.000109;
1142  YNOeq[16][1] = 0.001144;
1143  KdYNO[16][1] = 0.000754;
1144  YNOeq[16][2] = 0.001864;
1145  KdYNO[16][2] = 0.001323;
1146  YNOeq[16][3] = 0.002299;
1147  KdYNO[16][3] = 0.001679;
1148  YNOeq[16][4] = 0.002630;
1149  KdYNO[16][4] = 0.001974;
1150  YNOeq[16][5] = 0.002961;
1151  KdYNO[16][5] = 0.002269;
1152  YNOeq[16][6] = 0.003291;
1153  KdYNO[16][6] = 0.002565;
1154  YNOeq[16][7] = 0.003658;
1155  KdYNO[16][7] = 0.003131;
1156  YNOeq[16][8] = 0.004045;
1157  KdYNO[16][8] = 0.003843;
1158  YNOeq[16][9] = 0.004328;
1159  KdYNO[16][9] = 0.004167;
1160  YNOeq[16][10] = 0.004590;
1161  KdYNO[16][10] = 0.004491;
1162  YNOeq[16][11] = 0.004820;
1163  KdYNO[16][11] = 0.004800;
1164  YNOeq[17][0] = 0.000305;
1165  KdYNO[17][0] = 0.000389;
1166  YNOeq[17][1] = 0.000871;
1167  KdYNO[17][1] = 0.001117;
1168  YNOeq[17][2] = 0.001512;
1169  KdYNO[17][2] = 0.001698;
1170  YNOeq[17][3] = 0.002012;
1171  KdYNO[17][3] = 0.002534;
1172  YNOeq[17][4] = 0.002495;
1173  KdYNO[17][4] = 0.003402;
1174  YNOeq[17][5] = 0.002978;
1175  KdYNO[17][5] = 0.004271;
1176  YNOeq[17][6] = 0.003460;
1177  KdYNO[17][6] = 0.005140;
1178  YNOeq[17][7] = 0.003943;
1179  KdYNO[17][7] = 0.006009;
1180  YNOeq[17][8] = 0.004425;
1181  KdYNO[17][8] = 0.006877;
1182  YNOeq[17][9] = 0.004795;
1183  KdYNO[17][9] = 0.007841;
1184  YNOeq[17][10] = 0.005123;
1185  KdYNO[17][10] = 0.008969;
1186  YNOeq[17][11] = 0.005441;
1187  KdYNO[17][11] = 0.010388;
1188  YNOeq[18][0] = 0.000340;
1189  KdYNO[18][0] = 0.000618;
1190  YNOeq[18][1] = 0.001487;
1191  KdYNO[18][1] = 0.003637;
1192  YNOeq[18][2] = 0.002386;
1193  KdYNO[18][2] = 0.006320;
1194  YNOeq[18][3] = 0.002942;
1195  KdYNO[18][3] = 0.007955;
1196  YNOeq[18][4] = 0.003358;
1197  KdYNO[18][4] = 0.009197;
1198  YNOeq[18][5] = 0.003774;
1199  KdYNO[18][5] = 0.010439;
1200  YNOeq[18][6] = 0.004190;
1201  KdYNO[18][6] = 0.011681;
1202  YNOeq[18][7] = 0.004658;
1203  KdYNO[18][7] = 0.014287;
1204  YNOeq[18][8] = 0.005154;
1205  KdYNO[18][8] = 0.017620;
1206  YNOeq[18][9] = 0.005490;
1207  KdYNO[18][9] = 0.018637;
1208  YNOeq[18][10] = 0.005836;
1209  KdYNO[18][10] = 0.020254;
1210  YNOeq[18][11] = 0.006110;
1211  KdYNO[18][11] = 0.021109;
1212  YNOeq[19][0] = 0.000498;
1213  KdYNO[19][0] = 0.002188;
1214  YNOeq[19][1] = 0.001192;
1215  KdYNO[19][1] = 0.005203;
1216  YNOeq[19][2] = 0.001950;
1217  KdYNO[19][2] = 0.007094;
1218  YNOeq[19][3] = 0.002567;
1219  KdYNO[19][3] = 0.010552;
1220  YNOeq[19][4] = 0.003165;
1221  KdYNO[19][4] = 0.014212;
1222  YNOeq[19][5] = 0.003763;
1223  KdYNO[19][5] = 0.017873;
1224  YNOeq[19][6] = 0.004361;
1225  KdYNO[19][6] = 0.021533;
1226  YNOeq[19][7] = 0.004959;
1227  KdYNO[19][7] = 0.025193;
1228  YNOeq[19][8] = 0.005557;
1229  KdYNO[19][8] = 0.028854;
1230  YNOeq[19][9] = 0.006029;
1231  KdYNO[19][9] = 0.033036;
1232  YNOeq[19][10] = 0.006430;
1233  KdYNO[19][10] = 0.037389;
1234  YNOeq[19][11] = 0.006826;
1235  KdYNO[19][11] = 0.042859;
1236  YNOeq[20][0] = 0.000532;
1237  KdYNO[20][0] = 0.003029;
1238  YNOeq[20][1] = 0.001892;
1239  KdYNO[20][1] = 0.015222;
1240  YNOeq[20][2] = 0.002987;
1241  KdYNO[20][2] = 0.026207;
1242  YNOeq[20][3] = 0.003679;
1243  KdYNO[20][3] = 0.032808;
1244  YNOeq[20][4] = 0.004194;
1245  KdYNO[20][4] = 0.037504;
1246  YNOeq[20][5] = 0.004708;
1247  KdYNO[20][5] = 0.042201;
1248  YNOeq[20][6] = 0.005223;
1249  KdYNO[20][6] = 0.046897;
1250  YNOeq[20][7] = 0.005806;
1251  KdYNO[20][7] = 0.057298;
1252  YNOeq[20][8] = 0.006425;
1253  KdYNO[20][8] = 0.070743;
1254  YNOeq[20][9] = 0.006822;
1255  KdYNO[20][9] = 0.073583;
1256  YNOeq[20][10] = 0.007264;
1257  KdYNO[20][10] = 0.080335;
1258  YNOeq[20][11] = 0.007589;
1259  KdYNO[20][11] = 0.082249;
1260  YNOeq[21][0] = 0.000769;
1261  KdYNO[21][0] = 0.010316;
1262  YNOeq[21][1] = 0.001603;
1263  KdYNO[21][1] = 0.021430;
1264  YNOeq[21][2] = 0.002483;
1265  KdYNO[21][2] = 0.027015;
1266  YNOeq[21][3] = 0.003225;
1267  KdYNO[21][3] = 0.039702;
1268  YNOeq[21][4] = 0.003949;
1269  KdYNO[21][4] = 0.053304;
1270  YNOeq[21][5] = 0.004673;
1271  KdYNO[21][5] = 0.066906;
1272  YNOeq[21][6] = 0.005397;
1273  KdYNO[21][6] = 0.080508;
1274  YNOeq[21][7] = 0.006121;
1275  KdYNO[21][7] = 0.094110;
1276  YNOeq[21][8] = 0.006845;
1277  KdYNO[21][8] = 0.107712;
1278  YNOeq[21][9] = 0.007431;
1279  KdYNO[21][9] = 0.123569;
1280  YNOeq[21][10] = 0.007914;
1281  KdYNO[21][10] = 0.138603;
1282  YNOeq[21][11] = 0.008398;
1283  KdYNO[21][11] = 0.157431;
1284  YNOeq[22][0] = 0.000733;
1285  KdYNO[22][0] = 0.008818;
1286  YNOeq[22][1] = 0.002364;
1287  KdYNO[22][1] = 0.055774;
1288  YNOeq[22][2] = 0.003664;
1289  KdYNO[22][2] = 0.095179;
1290  YNOeq[22][3] = 0.004507;
1291  KdYNO[22][3] = 0.118901;
1292  YNOeq[22][4] = 0.005133;
1293  KdYNO[22][4] = 0.135043;
1294  YNOeq[22][5] = 0.005759;
1295  KdYNO[22][5] = 0.151184;
1296  YNOeq[22][6] = 0.006385;
1297  KdYNO[22][6] = 0.167326;
1298  YNOeq[22][7] = 0.007095;
1299  KdYNO[22][7] = 0.203906;
1300  YNOeq[22][8] = 0.007850;
1301  KdYNO[22][8] = 0.251396;
1302  YNOeq[22][9] = 0.008317;
1303  KdYNO[22][9] = 0.258692;
1304  YNOeq[22][10] = 0.008866;
1305  KdYNO[22][10] = 0.283291;
1306  YNOeq[22][11] = 0.009251;
1307  KdYNO[22][11] = 0.286651;
1308  YNOeq[23][0] = 0.001135;
1309  KdYNO[23][0] = 0.041490;
1310  YNOeq[23][1] = 0.002114;
1311  KdYNO[23][1] = 0.078217;
1312  YNOeq[23][2] = 0.003119;
1313  KdYNO[23][2] = 0.093379;
1314  YNOeq[23][3] = 0.003993;
1315  KdYNO[23][3] = 0.135042;
1316  YNOeq[23][4] = 0.004850;
1317  KdYNO[23][4] = 0.180119;
1318  YNOeq[23][5] = 0.005708;
1319  KdYNO[23][5] = 0.225196;
1320  YNOeq[23][6] = 0.006566;
1321  KdYNO[23][6] = 0.270273;
1322  YNOeq[23][7] = 0.007423;
1323  KdYNO[23][7] = 0.315350;
1324  YNOeq[23][8] = 0.008281;
1325  KdYNO[23][8] = 0.360426;
1326  YNOeq[23][9] = 0.008991;
1327  KdYNO[23][9] = 0.413773;
1328  YNOeq[23][10] = 0.009566;
1329  KdYNO[23][10] = 0.460907;
1330  YNOeq[23][11] = 0.010147;
1331  KdYNO[23][11] = 0.519839;
1332  YNOeq[24][0] = 0.001060;
1333  KdYNO[24][0] = 0.032900;
1334  YNOeq[24][1] = 0.002782;
1335  KdYNO[24][1] = 0.147795;
1336  YNOeq[24][2] = 0.003941;
1337  KdYNO[24][2] = 0.167231;
1338  YNOeq[24][3] = 0.005066;
1339  KdYNO[24][3] = 0.258906;
1340  YNOeq[24][4] = 0.005704;
1341  KdYNO[24][4] = 0.264210;
1342  YNOeq[24][5] = 0.006342;
1343  KdYNO[24][5] = 0.269515;
1344  YNOeq[24][6] = 0.007100;
1345  KdYNO[24][6] = 0.328396;
1346  YNOeq[24][7] = 0.007866;
1347  KdYNO[24][7] = 0.390753;
1348  YNOeq[24][8] = 0.008632;
1349  KdYNO[24][8] = 0.453110;
1350  YNOeq[24][9] = 0.009542;
1351  KdYNO[24][9] = 0.635565;
1352  YNOeq[24][10] = 0.010467;
1353  KdYNO[24][10] = 0.829667;
1354  YNOeq[24][11] = 0.011082;
1355  KdYNO[24][11] = 0.901448;
1356  YNOeq[25][0] = 0.001260;
1357  KdYNO[25][0] = 0.060880;
1358  YNOeq[25][1] = 0.003050;
1359  KdYNO[25][1] = 0.229015;
1360  YNOeq[25][2] = 0.004197;
1361  KdYNO[25][2] = 0.310817;
1362  YNOeq[25][3] = 0.005218;
1363  KdYNO[25][3] = 0.375920;
1364  YNOeq[25][4] = 0.006240;
1365  KdYNO[25][4] = 0.441022;
1366  YNOeq[25][5] = 0.007151;
1367  KdYNO[25][5] = 0.533430;
1368  YNOeq[25][6] = 0.008052;
1369  KdYNO[25][6] = 0.717659;
1370  YNOeq[25][7] = 0.008954;
1371  KdYNO[25][7] = 0.901887;
1372  YNOeq[25][8] = 0.009856;
1373  KdYNO[25][8] = 1.086116;
1374  YNOeq[25][9] = 0.010694;
1375  KdYNO[25][9] = 1.250362;
1376  YNOeq[25][10] = 0.011370;
1377  KdYNO[25][10] = 1.385262;
1378  YNOeq[25][11] = 0.012055;
1379  KdYNO[25][11] = 1.553527;
1380  YNOeq[26][0] = 0.001407;
1381  KdYNO[26][0] = 0.092724;
1382  YNOeq[26][1] = 0.002875;
1383  KdYNO[26][1] = 0.212964;
1384  YNOeq[26][2] = 0.004343;
1385  KdYNO[26][2] = 0.333203;
1386  YNOeq[26][3] = 0.005701;
1387  KdYNO[26][3] = 0.555891;
1388  YNOeq[26][4] = 0.006614;
1389  KdYNO[26][4] = 0.700366;
1390  YNOeq[26][5] = 0.007426;
1391  KdYNO[26][5] = 0.738973;
1392  YNOeq[26][6] = 0.008332;
1393  KdYNO[26][6] = 0.913561;
1394  YNOeq[26][7] = 0.009244;
1395  KdYNO[26][7] = 1.096972;
1396  YNOeq[26][8] = 0.010156;
1397  KdYNO[26][8] = 1.280383;
1398  YNOeq[26][9] = 0.011238;
1399  KdYNO[26][9] = 1.810014;
1400  YNOeq[26][10] = 0.012336;
1401  KdYNO[26][10] = 2.373222;
1402  YNOeq[26][11] = 0.013061;
1403  KdYNO[26][11] = 2.573348;
1404  YNOeq[27][0] = 0.001732;
1405  KdYNO[27][0] = 0.191921;
1406  YNOeq[27][1] = 0.003699;
1407  KdYNO[27][1] = 0.638817;
1408  YNOeq[27][2] = 0.004995;
1409  KdYNO[27][2] = 0.858891;
1410  YNOeq[27][3] = 0.006161;
1411  KdYNO[27][3] = 1.035088;
1412  YNOeq[27][4] = 0.007328;
1413  KdYNO[27][4] = 1.211285;
1414  YNOeq[27][5] = 0.008372;
1415  KdYNO[27][5] = 1.458364;
1416  YNOeq[27][6] = 0.009428;
1417  KdYNO[27][6] = 1.965043;
1418  YNOeq[27][7] = 0.010484;
1419  KdYNO[27][7] = 2.471722;
1420  YNOeq[27][8] = 0.011540;
1421  KdYNO[27][8] = 2.978402;
1422  YNOeq[27][9] = 0.012519;
1423  KdYNO[27][9] = 3.428585;
1424  YNOeq[27][10] = 0.013015;
1425  KdYNO[27][10] = 3.273183;
1426  YNOeq[27][11] = 0.014101;
1427  KdYNO[27][11] = 4.229005;
1428  YNOeq[28][0] = 0.001892;
1429  KdYNO[28][0] = 0.271459;
1430  YNOeq[28][1] = 0.003479;
1431  KdYNO[28][1] = 0.548492;
1432  YNOeq[28][2] = 0.005065;
1433  KdYNO[28][2] = 0.825524;
1434  YNOeq[28][3] = 0.006468;
1435  KdYNO[28][3] = 1.347835;
1436  YNOeq[28][4] = 0.007661;
1437  KdYNO[28][4] = 1.792081;
1438  YNOeq[28][5] = 0.008605;
1439  KdYNO[28][5] = 1.891306;
1440  YNOeq[28][6] = 0.009663;
1441  KdYNO[28][6] = 2.357629;
1442  YNOeq[28][7] = 0.010729;
1443  KdYNO[28][7] = 2.847770;
1444  YNOeq[28][8] = 0.011795;
1445  KdYNO[28][8] = 3.337912;
1446  YNOeq[28][9] = 0.012961;
1447  KdYNO[28][9] = 4.311444;
1448  YNOeq[28][10] = 0.014040;
1449  KdYNO[28][10] = 5.344972;
1450  YNOeq[28][11] = 0.015074;
1451  KdYNO[28][11] = 6.473152;
1452  YNOeq[29][0] = 0.002311;
1453  KdYNO[29][0] = 0.547036;
1454  YNOeq[29][1] = 0.004436;
1455  KdYNO[29][1] = 1.623547;
1456  YNOeq[29][2] = 0.005878;
1457  KdYNO[29][2] = 2.167642;
1458  YNOeq[29][3] = 0.007187;
1459  KdYNO[29][3] = 2.608745;
1460  YNOeq[29][4] = 0.008496;
1461  KdYNO[29][4] = 3.049848;
1462  YNOeq[29][5] = 0.009585;
1463  KdYNO[29][5] = 3.545103;
1464  YNOeq[29][6] = 0.010636;
1465  KdYNO[29][6] = 4.567760;
1466  YNOeq[29][7] = 0.011811;
1467  KdYNO[29][7] = 5.626905;
1468  YNOeq[29][8] = 0.012986;
1469  KdYNO[29][8] = 6.686050;
1470  YNOeq[29][9] = 0.014160;
1471  KdYNO[29][9] = 7.745196;
1472  YNOeq[29][10] = 0.014952;
1473  KdYNO[29][10] = 7.970170;
1474  YNOeq[29][11] = 0.015965;
1475  KdYNO[29][11] = 9.430666;
1476  YNOeq[30][0] = 0.002479;
1477  KdYNO[30][0] = 0.729234;
1478  YNOeq[30][1] = 0.004172;
1479  KdYNO[30][1] = 1.328519;
1480  YNOeq[30][2] = 0.005864;
1481  KdYNO[30][2] = 1.927803;
1482  YNOeq[30][3] = 0.007296;
1483  KdYNO[30][3] = 2.736098;
1484  YNOeq[30][4] = 0.008593;
1485  KdYNO[30][4] = 3.615589;
1486  YNOeq[30][5] = 0.009854;
1487  KdYNO[30][5] = 4.471673;
1488  YNOeq[30][6] = 0.011080;
1489  KdYNO[30][6] = 5.647474;
1490  YNOeq[30][7] = 0.012303;
1491  KdYNO[30][7] = 6.844018;
1492  YNOeq[30][8] = 0.013526;
1493  KdYNO[30][8] = 8.040562;
1494  YNOeq[30][9] = 0.014740;
1495  KdYNO[30][9] = 9.826371;
1496  YNOeq[30][10] = 0.016016;
1497  KdYNO[30][10] = 12.394166;
1498  YNOeq[30][11] = 0.015965;
1499  KdYNO[30][11] = 9.430666;
1500  YNOeq[31][0] = 0.002998;
1501  KdYNO[31][0] = 1.423809;
1502  YNOeq[31][1] = 0.004828;
1503  KdYNO[31][1] = 2.519829;
1504  YNOeq[31][2] = 0.006421;
1505  KdYNO[31][2] = 3.574074;
1506  YNOeq[31][3] = 0.007964;
1507  KdYNO[31][3] = 4.619481;
1508  YNOeq[31][4] = 0.009551;
1509  KdYNO[31][4] = 6.218774;
1510  YNOeq[31][5] = 0.011081;
1511  KdYNO[31][5] = 8.172902;
1512  YNOeq[31][6] = 0.012270;
1513  KdYNO[31][6] = 10.093628;
1514  YNOeq[31][7] = 0.013460;
1515  KdYNO[31][7] = 12.014354;
1516  YNOeq[31][8] = 0.014649;
1517  KdYNO[31][8] = 13.935080;
1518  YNOeq[31][9] = 0.015837;
1519  KdYNO[31][9] = 15.873938;
1520  YNOeq[31][10] = 0.016990;
1521  KdYNO[31][10] = 18.181945;
1522  YNOeq[31][11] = 0.015965;
1523  KdYNO[31][11] = 9.430666;
1524  YNOeq[32][0] = 0.002987;
1525  KdYNO[32][0] = 1.478557;
1526  YNOeq[32][1] = 0.004827;
1527  KdYNO[32][1] = 2.797413;
1528  YNOeq[32][2] = 0.006595;
1529  KdYNO[32][2] = 4.249854;
1530  YNOeq[32][3] = 0.008122;
1531  KdYNO[32][3] = 6.151483;
1532  YNOeq[32][4] = 0.009650;
1533  KdYNO[32][4] = 8.053111;
1534  YNOeq[32][5] = 0.011178;
1535  KdYNO[32][5] = 9.954740;
1536  YNOeq[32][6] = 0.012568;
1537  KdYNO[32][6] = 12.597257;
1538  YNOeq[32][7] = 0.013949;
1539  KdYNO[32][7] = 15.287844;
1540  YNOeq[32][8] = 0.015331;
1541  KdYNO[32][8] = 17.978430;
1542  YNOeq[32][9] = 0.016690;
1543  KdYNO[32][9] = 21.877635;
1544  YNOeq[32][10] = 0.016990;
1545  KdYNO[32][10] = 18.181945;
1546  YNOeq[32][11] = 0.015965;
1547  KdYNO[32][11] = 9.430666;
1548  YNOeq[33][0] = 0.003794;
1549  KdYNO[33][0] = 3.423591;
1550  YNOeq[33][1] = 0.005713;
1551  KdYNO[33][1] = 5.641213;
1552  YNOeq[33][2] = 0.007425;
1553  KdYNO[33][2] = 7.845004;
1554  YNOeq[33][3] = 0.009093;
1555  KdYNO[33][3] = 10.045869;
1556  YNOeq[33][4] = 0.010826;
1557  KdYNO[33][4] = 13.462808;
1558  YNOeq[33][5] = 0.012517;
1559  KdYNO[33][5] = 17.651895;
1560  YNOeq[33][6] = 0.013877;
1561  KdYNO[33][6] = 21.407047;
1562  YNOeq[33][7] = 0.015237;
1563  KdYNO[33][7] = 25.162198;
1564  YNOeq[33][8] = 0.016553;
1565  KdYNO[33][8] = 28.970749;
1566  YNOeq[33][9] = 0.017734;
1567  KdYNO[33][9] = 32.850350;
1568  YNOeq[33][10] = 0.016990;
1569  KdYNO[33][10] = 18.181945;
1570  YNOeq[33][11] = 0.015965;
1571  KdYNO[33][11] = 9.430666;
1572  YNOeq[34][0] = 0.003697;
1573  KdYNO[34][0] = 3.302310;
1574  YNOeq[34][1] = 0.005641;
1575  KdYNO[34][1] = 5.913553;
1576  YNOeq[34][2] = 0.007524;
1577  KdYNO[34][2] = 8.840037;
1578  YNOeq[34][3] = 0.009206;
1579  KdYNO[34][3] = 12.826540;
1580  YNOeq[34][4] = 0.010888;
1581  KdYNO[34][4] = 16.813044;
1582  YNOeq[34][5] = 0.012570;
1583  KdYNO[34][5] = 20.799547;
1584  YNOeq[34][6] = 0.014116;
1585  KdYNO[34][6] = 26.307950;
1586  YNOeq[34][7] = 0.015652;
1587  KdYNO[34][7] = 31.915095;
1588  YNOeq[34][8] = 0.017189;
1589  KdYNO[34][8] = 37.522241;
1590  YNOeq[34][9] = 0.017734;
1591  KdYNO[34][9] = 32.850350;
1592  YNOeq[34][10] = 0.016990;
1593  KdYNO[34][10] = 18.181945;
1594  YNOeq[34][11] = 0.015965;
1595  KdYNO[34][11] = 9.430666;
1596  YNOeq[35][0] = 0.004045;
1597  KdYNO[35][0] = 4.577814;
1598  YNOeq[35][1] = 0.006250;
1599  KdYNO[35][1] = 9.798461;
1600  YNOeq[35][2] = 0.008273;
1601  KdYNO[35][2] = 14.855960;
1602  YNOeq[35][3] = 0.010156;
1603  KdYNO[35][3] = 19.786978;
1604  YNOeq[35][4] = 0.012162;
1605  KdYNO[35][4] = 26.201035;
1606  YNOeq[35][5] = 0.013852;
1607  KdYNO[35][5] = 33.581057;
1608  YNOeq[35][6] = 0.015309;
1609  KdYNO[35][6] = 39.701992;
1610  YNOeq[35][7] = 0.016871;
1611  KdYNO[35][7] = 48.099433;
1612  YNOeq[35][8] = 0.018213;
1613  KdYNO[35][8] = 54.505888;
1614  YNOeq[35][9] = 0.017734;
1615  KdYNO[35][9] = 32.850350;
1616  YNOeq[35][10] = 0.016990;
1617  KdYNO[35][10] = 18.181945;
1618  YNOeq[35][11] = 0.015965;
1619  KdYNO[35][11] = 9.430666;
1620  YNOeq[36][0] = 0.004503;
1621  KdYNO[36][0] = 7.015211;
1622  YNOeq[36][1] = 0.006540;
1623  KdYNO[36][1] = 11.921544;
1624  YNOeq[36][2] = 0.008530;
1625  KdYNO[36][2] = 17.499392;
1626  YNOeq[36][3] = 0.010358;
1627  KdYNO[36][3] = 25.335258;
1628  YNOeq[36][4] = 0.012187;
1629  KdYNO[36][4] = 33.171123;
1630  YNOeq[36][5] = 0.014016;
1631  KdYNO[36][5] = 41.006989;
1632  YNOeq[36][6] = 0.015687;
1633  KdYNO[36][6] = 50.802898;
1634  YNOeq[36][7] = 0.017470;
1635  KdYNO[36][7] = 59.199654;
1636  YNOeq[36][8] = 0.018213;
1637  KdYNO[36][8] = 54.505888;
1638  YNOeq[36][9] = 0.017734;
1639  KdYNO[36][9] = 32.850350;
1640  YNOeq[36][10] = 0.016990;
1641  KdYNO[36][10] = 18.181945;
1642  YNOeq[36][11] = 0.015965;
1643  KdYNO[36][11] = 9.430666;
1644  YNOeq[37][0] = 0.004910;
1645  KdYNO[37][0] = 9.582819;
1646  YNOeq[37][1] = 0.007222;
1647  KdYNO[37][1] = 19.290391;
1648  YNOeq[37][2] = 0.009365;
1649  KdYNO[37][2] = 28.791551;
1650  YNOeq[37][3] = 0.011378;
1651  KdYNO[37][3] = 38.132688;
1652  YNOeq[37][4] = 0.013425;
1653  KdYNO[37][4] = 47.851436;
1654  YNOeq[37][5] = 0.015123;
1655  KdYNO[37][5] = 58.733176;
1656  YNOeq[37][6] = 0.016886;
1657  KdYNO[37][6] = 73.864675;
1658  YNOeq[37][7] = 0.018403;
1659  KdYNO[37][7] = 84.540857;
1660  YNOeq[37][8] = 0.018213;
1661  KdYNO[37][8] = 54.505888;
1662  YNOeq[37][9] = 0.017734;
1663  KdYNO[37][9] = 32.850350;
1664  YNOeq[37][10] = 0.016990;
1665  KdYNO[37][10] = 18.181945;
1666  YNOeq[37][11] = 0.015965;
1667  KdYNO[37][11] = 9.430666;
1668  YNOeq[38][0] = 0.005081;
1669  KdYNO[38][0] = 10.810942;
1670  YNOeq[38][1] = 0.007358;
1671  KdYNO[38][1] = 21.064593;
1672  YNOeq[38][2] = 0.009609;
1673  KdYNO[38][2] = 33.021948;
1674  YNOeq[38][3] = 0.011576;
1675  KdYNO[38][3] = 47.544379;
1676  YNOeq[38][4] = 0.013542;
1677  KdYNO[38][4] = 62.066810;
1678  YNOeq[38][5] = 0.015508;
1679  KdYNO[38][5] = 76.589241;
1680  YNOeq[38][6] = 0.017375;
1681  KdYNO[38][6] = 92.178864;
1682  YNOeq[38][7] = 0.018403;
1683  KdYNO[38][7] = 84.540857;
1684  YNOeq[38][8] = 0.018213;
1685  KdYNO[38][8] = 54.505888;
1686  YNOeq[38][9] = 0.017734;
1687  KdYNO[38][9] = 32.850350;
1688  YNOeq[38][10] = 0.016990;
1689  KdYNO[38][10] = 18.181945;
1690  YNOeq[38][11] = 0.015965;
1691  KdYNO[38][11] = 9.430666;
1692  YNOeq[39][0] = 0.005869;
1693  KdYNO[39][0] = 18.928014;
1694  YNOeq[39][1] = 0.008278;
1695  KdYNO[39][1] = 36.111013;
1696  YNOeq[39][2] = 0.010530;
1697  KdYNO[39][2] = 53.043002;
1698  YNOeq[39][3] = 0.012659;
1699  KdYNO[39][3] = 69.780392;
1700  YNOeq[39][4] = 0.014779;
1701  KdYNO[39][4] = 85.511999;
1702  YNOeq[39][5] = 0.016443;
1703  KdYNO[39][5] = 98.841862;
1704  YNOeq[39][6] = 0.018292;
1705  KdYNO[39][6] = 122.897344;
1706  YNOeq[39][7] = 0.018403;
1707  KdYNO[39][7] = 84.540857;
1708  YNOeq[39][8] = 0.018213;
1709  KdYNO[39][8] = 54.505888;
1710  YNOeq[39][9] = 0.017734;
1711  KdYNO[39][9] = 32.850350;
1712  YNOeq[39][10] = 0.016990;
1713  KdYNO[39][10] = 18.181945;
1714  YNOeq[39][11] = 0.015965;
1715  KdYNO[39][11] = 9.430666;
1716  YNOeq[40][0] = 0.006016;
1717  KdYNO[40][0] = 20.651069;
1718  YNOeq[40][1] = 0.008391;
1719  KdYNO[40][1] = 38.537591;
1720  YNOeq[40][2] = 0.010860;
1721  KdYNO[40][2] = 58.400257;
1722  YNOeq[40][3] = 0.013135;
1723  KdYNO[40][3] = 82.521261;
1724  YNOeq[40][4] = 0.015215;
1725  KdYNO[40][4] = 111.901382;
1726  YNOeq[40][5] = 0.017126;
1727  KdYNO[40][5] = 131.273590;
1728  YNOeq[40][6] = 0.018292;
1729  KdYNO[40][6] = 122.897344;
1730  YNOeq[40][7] = 0.018403;
1731  KdYNO[40][7] = 84.540857;
1732  YNOeq[40][8] = 0.018213;
1733  KdYNO[40][8] = 54.505888;
1734  YNOeq[40][9] = 0.017734;
1735  KdYNO[40][9] = 32.850350;
1736  YNOeq[40][10] = 0.016990;
1737  KdYNO[40][10] = 18.181945;
1738  YNOeq[40][11] = 0.015965;
1739  KdYNO[40][11] = 9.430666;
1740  YNOeq[41][0] = 0.006678;
1741  KdYNO[41][0] = 30.590106;
1742  YNOeq[41][1] = 0.009003;
1743  KdYNO[41][1] = 56.851468;
1744  YNOeq[41][2] = 0.011337;
1745  KdYNO[41][2] = 85.428380;
1746  YNOeq[41][3] = 0.013641;
1747  KdYNO[41][3] = 112.821628;
1748  YNOeq[41][4] = 0.015914;
1749  KdYNO[41][4] = 139.021404;
1750  YNOeq[41][5] = 0.017822;
1751  KdYNO[41][5] = 162.375287;
1752  YNOeq[41][6] = 0.018292;
1753  KdYNO[41][6] = 122.897344;
1754  YNOeq[41][7] = 0.018403;
1755  KdYNO[41][7] = 84.540857;
1756  YNOeq[41][8] = 0.018213;
1757  KdYNO[41][8] = 54.505888;
1758  YNOeq[41][9] = 0.017734;
1759  KdYNO[41][9] = 32.850350;
1760  YNOeq[41][10] = 0.016990;
1761  KdYNO[41][10] = 18.181945;
1762  YNOeq[41][11] = 0.015965;
1763  KdYNO[41][11] = 9.430666;
1764  YNOeq[42][0] = 0.007042;
1765  KdYNO[42][0] = 37.697174;
1766  YNOeq[42][1] = 0.009506;
1767  KdYNO[42][1] = 67.480598;
1768  YNOeq[42][2] = 0.012090;
1769  KdYNO[42][2] = 100.891368;
1770  YNOeq[42][3] = 0.014310;
1771  KdYNO[42][3] = 126.744041;
1772  YNOeq[42][4] = 0.016455;
1773  KdYNO[42][4] = 166.893523;
1774  YNOeq[42][5] = 0.017822;
1775  KdYNO[42][5] = 162.375287;
1776  YNOeq[42][6] = 0.018292;
1777  KdYNO[42][6] = 122.897344;
1778  YNOeq[42][7] = 0.018403;
1779  KdYNO[42][7] = 84.540857;
1780  YNOeq[42][8] = 0.018213;
1781  KdYNO[42][8] = 54.505888;
1782  YNOeq[42][9] = 0.017734;
1783  KdYNO[42][9] = 32.850350;
1784  YNOeq[42][10] = 0.016990;
1785  KdYNO[42][10] = 18.181945;
1786  YNOeq[42][11] = 0.015965;
1787  KdYNO[42][11] = 9.430666;
1788  YNOeq[43][0] = 0.007797;
1789  KdYNO[43][0] = 54.697991;
1790  YNOeq[43][1] = 0.010132;
1791  KdYNO[43][1] = 94.655766;
1792  YNOeq[43][2] = 0.012471;
1793  KdYNO[43][2] = 137.909658;
1794  YNOeq[43][3] = 0.014842;
1795  KdYNO[43][3] = 175.933649;
1796  YNOeq[43][4] = 0.017201;
1797  KdYNO[43][4] = 214.977433;
1798  YNOeq[43][5] = 0.017822;
1799  KdYNO[43][5] = 162.375287;
1800  YNOeq[43][6] = 0.018292;
1801  KdYNO[43][6] = 122.897344;
1802  YNOeq[43][7] = 0.018403;
1803  KdYNO[43][7] = 84.540857;
1804  YNOeq[43][8] = 0.018213;
1805  KdYNO[43][8] = 54.505888;
1806  YNOeq[43][9] = 0.017734;
1807  KdYNO[43][9] = 32.850350;
1808  YNOeq[43][10] = 0.016990;
1809  KdYNO[43][10] = 18.181945;
1810  YNOeq[43][11] = 0.015965;
1811  KdYNO[43][11] = 9.430666;
1812  YNOeq[44][0] = 0.008082;
1813  KdYNO[44][0] = 62.810606;
1814  YNOeq[44][1] = 0.010710;
1815  KdYNO[44][1] = 114.217438;
1816  YNOeq[44][2] = 0.013385;
1817  KdYNO[44][2] = 168.602766;
1818  YNOeq[44][3] = 0.015685;
1819  KdYNO[44][3] = 210.612953;
1820  YNOeq[44][4] = 0.017201;
1821  KdYNO[44][4] = 214.977433;
1822  YNOeq[44][5] = 0.017822;
1823  KdYNO[44][5] = 162.375287;
1824  YNOeq[44][6] = 0.018292;
1825  KdYNO[44][6] = 122.897344;
1826  YNOeq[44][7] = 0.018403;
1827  KdYNO[44][7] = 84.540857;
1828  YNOeq[44][8] = 0.018213;
1829  KdYNO[44][8] = 54.505888;
1830  YNOeq[44][9] = 0.017734;
1831  KdYNO[44][9] = 32.850350;
1832  YNOeq[44][10] = 0.016990;
1833  KdYNO[44][10] = 18.181945;
1834  YNOeq[44][11] = 0.015965;
1835  KdYNO[44][11] = 9.430666;
1836  YNOeq[45][0] = 0.008989;
1837  KdYNO[45][0] = 94.851480;
1838  YNOeq[45][1] = 0.011436;
1839  KdYNO[45][1] = 142.583015;
1840  YNOeq[45][2] = 0.013888;
1841  KdYNO[45][2] = 197.204589;
1842  YNOeq[45][3] = 0.016267;
1843  KdYNO[45][3] = 262.823904;
1844  YNOeq[45][4] = 0.017201;
1845  KdYNO[45][4] = 214.977433;
1846  YNOeq[45][5] = 0.017822;
1847  KdYNO[45][5] = 162.375287;
1848  YNOeq[45][6] = 0.018292;
1849  KdYNO[45][6] = 122.897344;
1850  YNOeq[45][7] = 0.018403;
1851  KdYNO[45][7] = 84.540857;
1852  YNOeq[45][8] = 0.018213;
1853  KdYNO[45][8] = 54.505888;
1854  YNOeq[45][9] = 0.017734;
1855  KdYNO[45][9] = 32.850350;
1856  YNOeq[45][10] = 0.016990;
1857  KdYNO[45][10] = 18.181945;
1858  YNOeq[45][11] = 0.015965;
1859  KdYNO[45][11] = 9.430666;
1860  YNOeq[46][0] = 0.009261;
1861  KdYNO[46][0] = 107.006637;
1862  YNOeq[46][1] = 0.011927;
1863  KdYNO[46][1] = 185.033724;
1864  YNOeq[46][2] = 0.014583;
1865  KdYNO[46][2] = 260.495804;
1866  YNOeq[46][3] = 0.016267;
1867  KdYNO[46][3] = 262.823904;
1868  YNOeq[46][4] = 0.017201;
1869  KdYNO[46][4] = 214.977433;
1870  YNOeq[46][5] = 0.017822;
1871  KdYNO[46][5] = 162.375287;
1872  YNOeq[46][6] = 0.018292;
1873  KdYNO[46][6] = 122.897344;
1874  YNOeq[46][7] = 0.018403;
1875  KdYNO[46][7] = 84.540857;
1876  YNOeq[46][8] = 0.018213;
1877  KdYNO[46][8] = 54.505888;
1878  YNOeq[46][9] = 0.017734;
1879  KdYNO[46][9] = 32.850350;
1880  YNOeq[46][10] = 0.016990;
1881  KdYNO[46][10] = 18.181945;
1882  YNOeq[46][11] = 0.015965;
1883  KdYNO[46][11] = 9.430666;
1884  YNOeq[47][0] = 0.010112;
1885  KdYNO[47][0] = 150.669136;
1886  YNOeq[47][1] = 0.012741;
1887  KdYNO[47][1] = 233.387507;
1888  YNOeq[47][2] = 0.015179;
1889  KdYNO[47][2] = 305.763808;
1890  YNOeq[47][3] = 0.016267;
1891  KdYNO[47][3] = 262.823904;
1892  YNOeq[47][4] = 0.017201;
1893  KdYNO[47][4] = 214.977433;
1894  YNOeq[47][5] = 0.017822;
1895  KdYNO[47][5] = 162.375287;
1896  YNOeq[47][6] = 0.018292;
1897  KdYNO[47][6] = 122.897344;
1898  YNOeq[47][7] = 0.018403;
1899  KdYNO[47][7] = 84.540857;
1900  YNOeq[47][8] = 0.018213;
1901  KdYNO[47][8] = 54.505888;
1902  YNOeq[47][9] = 0.017734;
1903  KdYNO[47][9] = 32.850350;
1904  YNOeq[47][10] = 0.016990;
1905  KdYNO[47][10] = 18.181945;
1906  YNOeq[47][11] = 0.015965;
1907  KdYNO[47][11] = 9.430666;
1908  YNOeq[48][0] = 0.010506;
1909  KdYNO[48][0] = 177.128602;
1910  YNOeq[48][1] = 0.013232;
1911  KdYNO[48][1] = 270.988925;
1912  YNOeq[48][2] = 0.015179;
1913  KdYNO[48][2] = 305.763808;
1914  YNOeq[48][3] = 0.016267;
1915  KdYNO[48][3] = 262.823904;
1916  YNOeq[48][4] = 0.017201;
1917  KdYNO[48][4] = 214.977433;
1918  YNOeq[48][5] = 0.017822;
1919  KdYNO[48][5] = 162.375287;
1920  YNOeq[48][6] = 0.018292;
1921  KdYNO[48][6] = 122.897344;
1922  YNOeq[48][7] = 0.018403;
1923  KdYNO[48][7] = 84.540857;
1924  YNOeq[48][8] = 0.018213;
1925  KdYNO[48][8] = 54.505888;
1926  YNOeq[48][9] = 0.017734;
1927  KdYNO[48][9] = 32.850350;
1928  YNOeq[48][10] = 0.016990;
1929  KdYNO[48][10] = 18.181945;
1930  YNOeq[48][11] = 0.015965;
1931  KdYNO[48][11] = 9.430666;
1932  YNOeq[49][0] = 0.011414;
1933  KdYNO[49][0] = 247.150173;
1934  YNOeq[49][1] = 0.013859;
1935  KdYNO[49][1] = 270.988925;
1936  YNOeq[49][2] = 0.015179;
1937  KdYNO[49][2] = 305.763808;
1938  YNOeq[49][3] = 0.016267;
1939  KdYNO[49][3] = 262.823904;
1940  YNOeq[49][4] = 0.017201;
1941  KdYNO[49][4] = 214.977433;
1942  YNOeq[49][5] = 0.017822;
1943  KdYNO[49][5] = 162.375287;
1944  YNOeq[49][6] = 0.018292;
1945  KdYNO[49][6] = 122.897344;
1946  YNOeq[49][7] = 0.018403;
1947  KdYNO[49][7] = 84.540857;
1948  YNOeq[49][8] = 0.018213;
1949  KdYNO[49][8] = 54.505888;
1950  YNOeq[49][9] = 0.017734;
1951  KdYNO[49][9] = 32.850350;
1952  YNOeq[49][10] = 0.016990;
1953  KdYNO[49][10] = 18.181945;
1954  YNOeq[49][11] = 0.015965;
1955  KdYNO[49][11] = 9.430666;
1956  YNOeq[50][0] = 0.011781;
1957  KdYNO[50][0] = 247.150173;
1958  YNOeq[50][1] = 0.013859;
1959  KdYNO[50][1] = 270.988925;
1960  YNOeq[50][2] = 0.015179;
1961  KdYNO[50][2] = 305.763808;
1962  YNOeq[50][3] = 0.016267;
1963  KdYNO[50][3] = 262.823904;
1964  YNOeq[50][4] = 0.017201;
1965  KdYNO[50][4] = 214.977433;
1966  YNOeq[50][5] = 0.017822;
1967  KdYNO[50][5] = 162.375287;
1968  YNOeq[50][6] = 0.018292;
1969  KdYNO[50][6] = 122.897344;
1970  YNOeq[50][7] = 0.018403;
1971  KdYNO[50][7] = 84.540857;
1972  YNOeq[50][8] = 0.018213;
1973  KdYNO[50][8] = 54.505888;
1974  YNOeq[50][9] = 0.017734;
1975  KdYNO[50][9] = 32.850350;
1976  YNOeq[50][10] = 0.016990;
1977  KdYNO[50][10] = 18.181945;
1978  YNOeq[50][11] = 0.015965;
1979  KdYNO[50][11] = 9.430666;
1980  YNOeq[51][0] = 0.012459;
1981  KdYNO[51][0] = 247.150173;
1982  YNOeq[51][1] = 0.013859;
1983  KdYNO[51][1] = 270.988925;
1984  YNOeq[51][2] = 0.015179;
1985  KdYNO[51][2] = 305.763808;
1986  YNOeq[51][3] = 0.016267;
1987  KdYNO[51][3] = 262.823904;
1988  YNOeq[51][4] = 0.017201;
1989  KdYNO[51][4] = 214.977433;
1990  YNOeq[51][5] = 0.017822;
1991  KdYNO[51][5] = 162.375287;
1992  YNOeq[51][6] = 0.018292;
1993  KdYNO[51][6] = 122.897344;
1994  YNOeq[51][7] = 0.018403;
1995  KdYNO[51][7] = 84.540857;
1996  YNOeq[51][8] = 0.018213;
1997  KdYNO[51][8] = 54.505888;
1998  YNOeq[51][9] = 0.017734;
1999  KdYNO[51][9] = 32.850350;
2000  YNOeq[51][10] = 0.016990;
2001  KdYNO[51][10] = 18.181945;
2002  YNOeq[51][11] = 0.015965;
2003  KdYNO[51][11] = 9.430666;
2004 
2005  if(temperature < 1250) {
2006  i = 1;
2007  } else if(temperature > 3750) {
2008  i = 50;
2009  } else {
2010  i = (temperature - 1200) / 50;
2011  }
2012  if(mO2 / mtotal <= 0.11) {
2013  j = (mO2 / mtotal + 0.01) / 0.01;
2014  di = (temperature - 1200) / 50 - i;
2015  dj = (mO2 / mtotal + 0.01) / 0.01 - j;
2016  i--;
2017  j--;
2018  } else {
2019  j = (0.10 + 0.01) / 0.01;
2020  di = (temperature - 1200) / 50 - i;
2021  dj = (0.10 + 0.01) / 0.01 - j;
2022  i--;
2023  j--;
2024  }
2025 
2026  if(i < 0) {
2027  i = 0;
2028  di = 0.00001;
2029  }
2030  if(j < 0) {
2031  j = j;
2032  dj = 0.00001;
2033  }
2034 
2035  if(di == 0)
2036  di = 0.00001;
2037  if(dj == 0)
2038  dj = 0.00001;
2039 
2040  *YNOeq_value = (1 / dj * (1 / di * YNOeq[i][j] + 1 / (1 - di) * YNOeq[i + 1][j]) / (1 / di + 1 / (1 - di)) + 1 /
2041  (1 - dj) * (1 / di * YNOeq[i][j + 1] + 1 / (1 - di) * YNOeq[i + 1][j + 1]) / (1 / di + 1 / (1 - di))) / (1 / dj + 1 /
2042  (1 - dj));
2043 
2044  *KdYNO_value = (1 / dj * (1 / di * KdYNO[i][j] + 1 / (1 - di) * KdYNO[i + 1][j]) / (1 / di + 1 / (1 - di)) + 1 /
2045  (1 - dj) * (1 / di * KdYNO[i][j + 1] + 1 / (1 - di) * KdYNO[i + 1][j + 1]) / (1 / di + 1 / (1 - di))) / (1 / dj + 1 /
2046  (1 - dj));
2047 
2048 }
2049 
2050 void FUNCTION_SOOT_C(double *soot_pre, double element_FI) {
2051 
2052  int i = 0;
2053  int di = 0, dk = 0;
2054  double *soot_pre_vector;
2055  double *FI_vector;
2056 
2057  soot_pre_vector = (double*) malloc(11 * sizeof(double));
2058  FI_vector = (double*) malloc(11 * sizeof(double));
2059 
2060  soot_pre_vector[0] = 0;
2061  soot_pre_vector[1] = 0.01;
2062  soot_pre_vector[2] = 0.02;
2063  soot_pre_vector[3] = 0.05;
2064  soot_pre_vector[4] = 0.1;
2065  soot_pre_vector[5] = 0.2;
2066  soot_pre_vector[6] = 0.6;
2067  soot_pre_vector[7] = 0.7;
2068  soot_pre_vector[8] = 0.8;
2069  soot_pre_vector[9] = 0.9;
2070  soot_pre_vector[10] = 1;
2071 
2072  FI_vector[0] = 0;
2073  FI_vector[1] = 1;
2074  FI_vector[2] = 1.5;
2075  FI_vector[3] = 2;
2076  FI_vector[4] = 2.2;
2077  FI_vector[5] = 2.5;
2078  FI_vector[6] = 3;
2079  FI_vector[7] = 3.5;
2080  FI_vector[8] = 4;
2081  FI_vector[9] = 5.5;
2082  FI_vector[10] = 10;
2083 
2084  if(element_FI >= 10) {
2085  *soot_pre = 1;
2086  } else {
2087  dk = 0;
2088  for(i = 0; i < 10; i++) {
2089  if((dk == 0) && (element_FI >= FI_vector[i])) {
2090  dk = 1;
2091  di = i;
2092  }
2093  }
2094  *soot_pre = (element_FI - FI_vector[di]) * (soot_pre_vector[di + 1] - soot_pre_vector[di]) /
2095  (FI_vector[di + 1] - FI_vector[di]) + soot_pre_vector[di];
2096 
2097  }
2098 
2099  free(soot_pre_vector);
2100  free(FI_vector);
2101 
2102 }
2103 
2104 // ****************************************************************************************************
2105 /* ______________________________________________________________________________________________________________________
2106 
2107  COMBUSTION MODEL MAIN FUNCTION
2108  ______________________________________________________________________________________________________________________ */
2109 
2110 // ****************************************************************************************************
2113 void ACT(double *engine_parameters, double *engine_model_constants, double *test_variables, double *injection_rate,
2114  double *CAD_injection_rate, int size_inlet_inj, int NIN, double *SOI, double *EOI, int CAI, double *CAD_exit,
2115  double *HRF_exit, double *ROHR_exit, double *p_cyl_exit, double *dp_da_cyl_exit, double *T_cyl_exit,
2116  double *H_cooler_exit, double *mean_var_exit, double *heat_transfer, double *injection_rate_exit,
2117  double *accum_injection_rate_exit, sINtype dataIN, sOUTtype *dataOUT) {
2118 
2119  // element mixture and combustion state
2120  enum stState {
2121  stLiquid = 0, stEvaporated = 1, stOvermixed = 2, stBurned_poor_premix = 3, stBurned_rich_premix = 4, stBurned_by_diffusion = 5, stBurned_by_second_diffusion = 6, stInto_diffusion_flame = 7
2122  };
2123 
2124  int elem_i = 0, pulso_i = 0;
2125 
2126  double realelementmfreac = 0.; // Auxiliar variable
2127 
2128  // General constants
2129  double PI = __cons::Pi;
2130  double Runiv = 8314.;
2131  FILE *fich, *foculto, *foculto2;
2132  double FRLOL = 0.; // Dosado en el lift-off del instante considerado
2133  double *Itot;
2134  int CalcRad[8] = { 0, 0, 0, 0, 0, 0, 0, 0 }; // Variable que indica si se ha de calcular o no la radiacion (al principio, NO).
2135  int RadCalc = 0; // Variable que indica el deseo del usuario de calcular o no la radiacion.
2136 
2137  double YO2aire = 0.23019;
2138  double YO2i = 0.; // Concentracion de O2 en el bowl en un instante
2139  double Zst = 0., Z1 = 0., Z0 = 0., TCC = 0., TSC = 0.; // Para calculo T (relaciones de estado)
2140  double mmH2O = 18.; // Molecular weight for H2O
2141  double mmCO2 = 44.; // Molecular weight for CO2
2142 
2143  char title[3000];
2144 
2145  // engine parameters assignation
2146 
2147  double Piston_D = engine_parameters[0]; /* Piston diamteter(m) */
2148  double S = engine_parameters[1]; /* Piston stroke(m) */
2149  double Crank_L = S / 2.; /* Crank length(m) */
2150  double Connecting_Rod_L = engine_parameters[3];
2151  /* Connecting rod length(m) */
2152  double E = engine_parameters[4]; /* Piston eccentricity(mm) */
2153  double Piston_Axis_D = engine_parameters[5]; /* Piston axis diameter(m) */
2154  double Piston_Crown_H = engine_parameters[6]; /* Piston crown height(m) */
2155  double DBowl = engine_parameters[7]; /* Maximum diameter of bowl(m) */
2156  double VBowl = engine_parameters[8]; /* Volume of Bowl(m3) */
2157  double M_Connecting_Rod = engine_parameters[9];
2158  /* Mass Connecting rod(Kg) */
2159  double M_P_R_PA = engine_parameters[10];
2160  /* Mass of piston+rings+piston axis(Kg) */
2161  double C_ESteel = engine_parameters[11];
2162  /* Elasticity module of the steel(N/m2) */
2163  double C_Mech_Defor = engine_parameters[12];
2164  /* Mechanical deformations coefficient */
2165  double C_MBLBY = engine_parameters[13];
2166  /* Coefficient of leak's section of Blow-by */
2167  double GCRatio = engine_parameters[14]; /* Geometric compresion ratio */
2168  double n_holes = engine_parameters[15]; /* Number of holes of the nozzle */
2169  double nozzle_d = engine_parameters[16]; /* nozzle diameter (m) */
2170  double dc = engine_parameters[17]; /* discharge coefficient of the nozzle */
2171  double CTM = engine_parameters[18]; /* parameter to calculate C1 */
2172  double WC1A = engine_parameters[19]; /* Constant A to calculate C1 */
2173  double WC1B = engine_parameters[20]; /* Constant B to calculate C1 */
2174  double C2 = engine_parameters[21];
2175  /* Coefficient to the Woschin's calculus */
2176  double IVC = engine_parameters[22]; /* Inlet valve closing (deg) */
2177  double EVO = engine_parameters[23]; /* Exhaust valve opening (deg) */
2178 
2179  double Kmixture1 = engine_parameters[24];
2180  /* mixture combustion model constants */
2181 
2182  // We check if the user wants to calculate Radiation
2183  if(Kmixture1 < 0.) {
2184  RadCalc = 1;
2185  Kmixture1 = -Kmixture1;
2186  } else
2187  RadCalc = 0;
2188 
2189  // combustion model constant assignation
2190 
2191  double Kmixture2 = engine_model_constants[0];
2192 
2193  double Kpmx1 = engine_model_constants[1];
2194  /* premix combustion model constants */
2195  double Kpmx2 = engine_model_constants[2];
2196  double Kpmx3 = engine_model_constants[3];
2197  double Kpmx4 = engine_model_constants[4];
2198  double Kpmx5 = engine_model_constants[5];
2199 
2200  double KID1 = engine_model_constants[6];
2201  /* Ignition delay time model constants */
2202  double KID2 = engine_model_constants[7];
2203  double KID3 = engine_model_constants[8];
2204  double KID4 = engine_model_constants[9];
2205  double KID5 = engine_model_constants[10];
2206 
2207  double KNOx1 = engine_model_constants[11];
2208  /* Nox emission prediction model constants */
2209  double KNOx2 = engine_model_constants[12];
2210 
2211  double EC = engine_model_constants[13];
2212 
2213  double KSOOTA1 = engine_model_constants[14];
2214  /* SOOT emission prediction model constants */
2215  double KSOOTA2 = engine_model_constants[15];
2216  double KSOOTA3 = engine_model_constants[16];
2217 
2218  double KSOOTA4 = engine_model_constants[17];
2219  double KSOOTA5 = engine_model_constants[18];
2220  double KSOOTA6 = engine_model_constants[19];
2221  double KSOOTA7 = engine_model_constants[20];
2222 
2223  double KSOOTC1 = dataIN.KSOOTC1;
2224  /* SOOT emission prediction model constants */
2225  // KSOOTC1=1.2771738e-3; // Imponemos este valor, para que todo cuadre...
2226  KSOOTC1 = 0.9e-3; // Imponemos este valor, para calar el soot...
2227 
2228  double mf_old = 0.; // Para poder calcular el dmf_reac de realelement.
2229 
2230  double SOOT_EVO_A = 0.;
2231  double SOOT_EVO_C = 0.;
2232 
2233  double Kswirl = 0; /* swirl correction model constant */
2234 
2235  // test variables assignement
2236 
2237  double speed = test_variables[0]; /* Engine speed (rpm) */
2238  double mairadm = test_variables[1];
2239  /* In-cylinder air mass at inlet valve closing (g) */
2240  double mairIVC = test_variables[2];
2241  /* In-cylinder air mass at inlet valve closing (g) */
2242  double temperatureIVC = test_variables[3];
2243  /* In_cylinder temperature at inlet valve closing(K) */
2244  double mfuel = test_variables[4]; /* Total injected fuel mass (mg) */
2245  double PCR = test_variables[5];
2246  double inlet_pres = test_variables[6]; /* Inlet pressure (Pa) */
2247  double exhaust_pres = test_variables[7]; /* Exhaust pressure (Pa) */
2248  double Cbb = test_variables[8]; /* Adjustment coefficient of Blow-by */
2249  double Atmosphere_press = test_variables[9]; /* Atmosphere pressure (Pa) */
2250  double inj_fuel_temp = test_variables[10];
2251  /* Injection fuel temperature (K) */
2252  double TCYL_HEAD = test_variables[11]; /* Piston head temperature(K) */
2253  double TCYL = test_variables[12]; /* Cylinder temperature(K) */
2254  double TPIS = test_variables[13]; /* Piston temperature(K) */
2255 
2256  double rateBG = 1. - mairadm / mairIVC;
2257 
2258  // fuel properties
2259 
2260  double rofuel = 830.; /* Fuel density (kg/m3) */
2261  double HC = 1.7843; /* Fuel H/C ratio */
2262  double AFe = 14.5; /* Stoichiometric air / fuel ratio */
2263  double mmfuel = 12. + HC; // Molecular weight of the fuel (for just 1 C atome)
2264  double Vc_factor = 0.99; /* Engine K factor */
2265  double D = 0.0018; /* Dissipation constant */
2266  double MW_air = 28.97; /* Molecular weight of air(g/mol) */
2267  double MW_fuel = 152.2; /* Molecular weight of fuel(g/mol) */
2268  double MW_burned = 29.13; /* Molecular weight of burned products(g/mol) */
2269  double HP = 42920000.; /* Net heat of combustion(J/Kg) */
2270  double T_Evaporation_fuel = 489.35;
2271  /* Temperature of evaporation of fuel(K) */
2272 
2273  // Dissipation constant is adapted depending on the engine speed
2274  if(speed > 2000.)
2275  D = D * speed / 2000.;
2276  // D=0.0018*2.5*speed/2000.; /* Dissipation constant */
2277 
2278  // individual variables declaration
2279 
2280  double PRECISION_ITERATION = 0.000005; /* Value to iterate */
2281  double delta_t = 0.; /* Calculus interval time (s) */
2282  int size = 0; /* Size of the temporal vector */
2283  double mf_burned = 0.; /* Fuel mass burned (-) */
2284  double mf_burned_pmx = 0.; // Fuel mass burned in premixed combustion.
2285  double Cylinder_capacity = 0.; /* Cylinder capacity(m3) */
2286  double VTDC = 0.; /* Volume at top dead center(m3) */
2287  double Gamma = 0.; /* Politropic exponent */
2288  double f = 0.; /* Real air fuel ratio */
2289  double mEGR = 0.; /* Burned gases mass(Kg) */
2290  double p_exit = 0.; /* In cylinder pressure at exhaust valve openning (bar) */
2291  double T_exit = 0.; /* In cylinder temperature at exhaust valve openning (K) */
2292  double WI_HP = 0; /* Indicated work of pressure high cicle(J) */
2293  double IMP_HP = 0; /* Indicated Mean Pressure of pressure high cicle(bar) */
2294  double pmax = 0; /* In-cylinder Maximum pressure(bar) */
2295  double Tmax = 0; /* In cylinder Maximum temperature(K) */
2296  double dp_da_max = 0; /* Maximum dp/d(alfa) (bar/deg) */
2297  double Uo_i = 0.; // Virtual velocity at instant i
2298  double YO2_bowl_i = 0.; // O2 mass fraction in the bowl at instant i
2299 
2300  // TEMPORAL VARIABLES DECLARATION
2301 
2302  double *acu_dmf; /* Accumulated injection rate (-) */
2303  double *time_vector; /* Time (s) [IVC - EVO] */
2304  double *time_vector_exit; /* Time (s) [-180 - +180] */
2305  double *CAD; /* Crank angle degrees (deg) [IVC - EVO] */
2306  double *complete_CAD; /* Crank angle degrees (deg) [-180 - +180] used to calculate
2307  the IMP_HP */
2308  double *roair; /* In-cylinder air density (kg/m3) */
2309  double *virtual_velocity; /* Virtual velocity (m/s) */
2310  double *inj_velocity; /* Injection velocity (m/s) */
2311  double *p_cyl; /* In-cylinder pressure (bar) [IVC - EVO] */
2312  double *dp_da_cyl;
2313  double *complete_p_cyl; /* In-cylinder pressure (bar) [-180 - +180] */
2314  double *T_cyl; /* Averaged in-cylinder temperature (K) */
2315  double *HRF; /* Heat release fraction (-) (FQL) */
2316  double *HRF_PMX; /* Heat release fraction for premixed combustion (-) */
2317  double *ROHR;
2318  double *Yair; /* Air fraction */
2319  double *Yfuel; /* Fuel fraction */
2320  double *Yburned; /* Stoichiometricly burned gas fraction */
2321  double *Mbb; /* Blow-by mass (Kg) */
2322  double *acu_Mbb; /* Accumulated blow-by mass (Kg) */
2323  double *U; /* Internal energy (J) */
2324  double *CV; /* At constant volume heat (J/Kg) */
2325  double *H_cooler; /* Heat hung over to the cooler (J/deg) */
2326  double *Qcylhead;
2327  double *Qcyl;
2328  double *Qpis;
2329  double *H; /* Convection coefficient */
2330  double *defor; /* Increase of volume due to cylinder pressure and inertia
2331  (m3) [IVC - EVO] */
2332  double *Rmixture; /* Constant of perfect gases */
2333  double *V_cyl; /* In-cylinder volume (m3) [IVC - EVO] */
2334  double *complete_V_cyl; /* In-cylinder volume (m3) [-180 - +180] used to shown the
2335  results */
2336  double *complete_deform;
2337  /* Increase of volume due to cylinder pressure and inertia (m3) [-180 - +180] used to calculate the IMP_HP */
2338  double *dmf;
2339  double *SOI_IM; /* start of injections of multiple pulse */
2340  double *EOI_IM; /* end of injections of multiple pulse */
2341  double *aux_vector;
2342  double *evol_mSoot; // Evolucion de la masa de soot en el chorro.
2343  double *mSootCil; // Evolucion de la masa de soot remanente en cilindro (viene del EGR).
2344  double *evol_Radiacion; // Evolucion del parametro de radiacion (m_soot*T^4)
2345 
2346  // ELEMENT AND SUB-ELEMENT VARIABLES DECLARATION
2347 
2348  int num_i = 0; /* Number of elements i */
2349  int num_j = 0; /* Number of sub-elements j */
2350  double *mixture_correction; /* Burnt limit for the sub-element j */
2351  int inj_num = 0;
2352  int inj_counter = 0;
2353  double rate_area = 0.;
2354 
2355  int i_aux = 0, j_aux = 0;
2356 
2357  // MULTIPLE INJECTION VARIABLES
2358 
2359  double **POI_IM;
2360  double **POC_IM;
2361  double **mfuel_i_IM;
2362  double **mfuel_ij_IM;
2363  int *num_i_IM;
2364  double *SOC_IM;
2365 
2366  // NOx emission variables
2367  double YNOeq_value = 0.;
2368  double KdYNO_value = 0.;
2369 
2370  double **YNOeq;
2371  double **KdYNO;
2372 
2373  // SOOT emission variables;
2374 
2375  double soot_pre = 0.;
2376  double RES_FSN = 0.;
2377 
2378  // ELEMENT TRACKING MATRIX: STATE, COMPOSITION AND CHARACTERISTICS
2379 
2380  struct stPropertiesElement {
2381  int state = 0; /* element mixture and combustion state */
2382  double mtotal = 0.; /* element total mass */
2383  double dmtotal = 0.; /* element mass include at mixture process */
2384  double mO2 = 0.; /* element O2 mass */
2385  double mN2 = 0.; /* element N2 mass */
2386  double mCO2 = 0.; /* element CO2 mass */
2387  double mH2O = 0.; /* element H2O mass */
2388  double mf_jet = 0.;
2389  /* element fuel mass injected according to injection rate law */
2390  double mf_reac = 0.; /* element fuel mass updated with mixture */
2391  double mf_evap = 0.;
2392  /* element fuel mass updated with mixture and combustion process */
2393  double C = 0.; /* element fuel concentration */
2394  double FI = 0.; /* element air-fuel ratio */
2395  double RID = 0.; /* element ignition delay time intensity */
2396  double Rpmx = 0.; /* element premix combustion intensity */
2397  double Rpmx_value = 0.;
2398  double mNOx = 0.; /* element NOx mass */
2399  double dNOx = 0.; /* element NOx mass increment */
2400  double HC = 0.; /* element HC mass */
2401  double CO = 0.; /* element CO mass */
2402  double mSOOT_A = 0.; /* element SOOT mass */
2403  double dSOOT_A = 0.;
2404  double mSOOT_B = 0.; /* element SOOT mass */
2405  double dSOOT_B = 0.;
2406  double mSOOT_C = 0.; /* element SOOT mass */
2407  double dSOOT_C = 0.;
2408  double TSD = 0.;
2409  double Tadib = 0.;
2410  double Pcyl_POC = 0.;
2411  double TNOx = 0.;
2412  double X = 0.;
2413  double soot_precursor = 0.;
2414  double tLOL = 0.; /* Time when the LOL is reached. */
2415  double FRLOL = 0.; /* Relative equiv. ratio at LOL. */
2416 
2417  }**element;
2418 
2419  double *XLO; // Longitud de lift-off
2420  double *PEN_min; // Punto mas retrasado del chorro
2421  double *PEN_max; // Punto mas adelantado del chorro
2422 
2423  struct stRealElementComposition {
2424  double mtotal = 0.; /* element total mass */
2425  double mO2 = 0.; /* element O2 mass */
2426  double mO = 0.; /* element O mass wich is provided by O2, H2O and CO2 */
2427  double mCH = 0.;
2428  double mN2 = 0.; /* element N2 mass */
2429  double mCO2 = 0.; /* element CO2 mass */
2430  double mH2O = 0.; /* element H2O mass */
2431  double mf_reac = 0.;
2432  /* element fuel mass updated with mixture and combustion process */
2433  double dmf_reac = 0.; /* change in element fuel mass at the time step */
2434  double mf_reac_pmx = 0.; /* element fuel mass for premixed combustion */
2435  double C = 0.; /* element fuel concentration */
2436  double FI = 0.; /* element air-fuel ratio */
2437  }**realelement;
2438 
2439  stControlElementComposition **elementcontrol;
2440 
2441  // inlet valve closing composition variables declaration
2442 
2443  double mtotal_IVC = 0.; /* Total mass at inlet valve closing */
2444  double mO2_IVC = 0.;
2445  /* O2 mass at inlet valve closing coming from air fresh mass */
2446  double mN2_IVC = 0.;
2447  /* N2 mass at inlet valve closing coming from air fresh mass */
2448  double mCO2_IVC = 0.;
2449  /* CO2 mass at inlet valve closing coming from exhaust gasses recirculated */
2450  double mH2O_IVC = 0.;
2451  /* H2O mass at inlet valve closing coming from exhaust gasses recirculated */
2452  double NOx_IVC = 0.;
2453  /* NOx mass at inlet valve closing coming from exhaust gasses recirculated emissions */
2454  double mSOOT_IVC_B = 0.;
2455  /* SOOT mass at inlet valve closing coming from exhaust gasses recirculated emissions */
2456  double mSOOT_IVC_A = 0.;
2457  double mSOOT_IVC_C = 0.;
2458  double CO_IVC = 0.;
2459  double HC_IVC = 0.;
2460  double YO2IVC = 0.; /* Mass oxygen concentration at inlet valve closing */
2461  double YN2IVC = 0.; /* Mass nitrogen concentration at inlet valve closing */
2462  double YCO2IVC = 0.;
2463  double YH2OIVC = 0.;
2464 
2465  // bowl variables definition: composition, species and characteristics
2466 
2467  double *mtotal_bowl;
2468  double *mO2_bowl;
2469  double *mN2_bowl;
2470  double *mCO2_bowl;
2471  double *mH2O_bowl;
2472  double *mHC_bowl;
2473  double *Tsq_cyl;
2474  double *Tadib_cyl;
2475  double *mNOx_bowl;
2476  double *HC_bowl;
2477  double *mCO_bowl;
2478  double *mSOOT_bowl_A;
2479  double *mSOOT_bowl_B;
2480  double *mSOOT_bowl_C;
2481 
2482  double mSOOT_bowl_A_i_burned = 0.;
2483  double mSOOT_bowl_B_i_burned = 0.;
2484  double mSOOT_bowl_C_i_burned = 0.;
2485 
2486  // Dead volume variables definition: composition, species and characteristics
2487 
2488  double *mtotal_Vc;
2489  double *mO2_Vc;
2490  double *mN2_Vc;
2491  double *mCO2_Vc;
2492  double *mH2O_Vc;
2493  double *mNOx_Vc;
2494  double *mSOOT_Vc_A;
2495  double *mSOOT_Vc_B;
2496  double *mSOOT_Vc_C;
2497  double *mHC_Vc;
2498  double *dmtotal_Gfactor;
2499  /* increment mass exchange between bowl and dead volume */
2500 
2501  // stoichometric constant link with HC fuel ratio
2502  double Kst1 = 0.;
2503  double Kst2 = 0.;
2504  double Kst3 = 0.;
2505  double Kst4 = 0.;
2506  double Kst5 = 0.;
2507  double Kst6 = 0.;
2508 
2511  int complete_size = 0;
2512  /* vector size from -180 to +180 with a delta_CAD increment */
2513  int complete_prev_size = 0;
2514  /* vector size from -180 to IVC with a delta_CAD increment */
2515  int complete_post_size = 0;
2516  /* vector size from EVO to +180 with a delta_CAD increment */
2517  double auxiliar = 0.;
2518  double *vector_to_interpolate;
2519 
2520  double counter_CAD_1 = 0.;
2521  double counter_CAD_2 = 0.;
2522  int counter = 0;
2523  int m = 0, i = 0, j = 0, aux = 0;
2524  double Y1 = 0., T1 = 0., RID1 = 0.;
2525 
2526  double maximum = 0.;
2527  double minimum = 0.;
2528  double aux_mfuel = 0., Tot = 0.;
2529 
2530  int element_value = 0;
2531 
2532  aux_mfuel = mfuel;
2533 
2534  if(mfuel == 0.) {
2535  mfuel = 1.e-6;
2536  aux_mfuel = 0.;
2537  }
2538 
2539  // crank angle degree and time vector creation: time vector is created with a constant increment 20 microseconds
2540  // vector size calculation with the name 'size'. this will be the dimension of all the instantaneous cylinder vectors
2541 
2546  double Ang_Grab = 0.;
2547  Ang_Grab = -180.0; // angulo en el que se grabaran los datos
2548 
2549  if(RadCalc == 1 && Ang_Grab > -179.) {
2550  foculto = fopen("elementos.csv", "w");
2551  if(foculto == NULL)
2552  exit(-1);
2553 
2554  foculto2 = fopen("data_eje.csv", "w");
2555  if(foculto2 == NULL)
2556  exit(-1);
2557  }
2558 
2561  delta_t = 2. * 1e-5;
2562  size = 0;
2563  auxiliar = IVC * 60. / (360. * speed);
2564  do {
2565  size = size + 1;
2566  auxiliar = auxiliar + delta_t;
2567  } while(auxiliar <= EVO * 60. / (360. * speed));
2568 
2569  CAD = (double*) malloc(size * sizeof(double));
2570  time_vector = (double*) malloc(size * sizeof(double));
2571 
2572  for(counter = 0; counter < size; counter++) {
2573  if(counter == 0) {
2574  time_vector[counter] = IVC * 60. / (360. * speed);
2575  CAD[counter] = IVC;
2576  } else {
2577  time_vector[counter] = time_vector[counter - 1] + delta_t;
2578  CAD[counter] = CAD[counter - 1] + delta_t * 360. * speed / 60.;
2579  }
2580  }
2581 
2582  // instantaneous cylinder vectors take format with size dimension
2583 
2584  SOI_IM = (double*) malloc(8 * sizeof(double));
2585  EOI_IM = (double*) malloc(8 * sizeof(double));
2586  SOC_IM = (double*) malloc(8 * sizeof(double));
2587 
2588  mtotal_bowl = (double*) malloc(size * sizeof(double));
2589  mO2_bowl = (double*) malloc(size * sizeof(double));
2590  mN2_bowl = (double*) malloc(size * sizeof(double));
2591  mCO2_bowl = (double*) malloc(size * sizeof(double));
2592  mH2O_bowl = (double*) malloc(size * sizeof(double));
2593  mHC_bowl = (double*) malloc(size * sizeof(double));
2594  Tsq_cyl = (double*) malloc(size * sizeof(double));
2595  Tadib_cyl = (double*) malloc(size * sizeof(double));
2596  mNOx_bowl = (double*) malloc(size * sizeof(double));
2597  mCO_bowl = (double*) malloc(size * sizeof(double));
2598  HC_bowl = (double*) malloc(size * sizeof(double));
2599  mSOOT_bowl_A = (double*) malloc(size * sizeof(double));
2600  mSOOT_bowl_B = (double*) malloc(size * sizeof(double));
2601  mSOOT_bowl_C = (double*) malloc(size * sizeof(double));
2602 
2603  mtotal_Vc = (double*) malloc(size * sizeof(double));
2604  mO2_Vc = (double*) malloc(size * sizeof(double));
2605  mN2_Vc = (double*) malloc(size * sizeof(double));
2606  mCO2_Vc = (double*) malloc(size * sizeof(double));
2607  mH2O_Vc = (double*) malloc(size * sizeof(double));
2608  mNOx_Vc = (double*) malloc(size * sizeof(double));
2609  mSOOT_Vc_A = (double*) malloc(size * sizeof(double));
2610  mSOOT_Vc_B = (double*) malloc(size * sizeof(double));
2611  mSOOT_Vc_C = (double*) malloc(size * sizeof(double));
2612  mHC_Vc = (double*) malloc(size * sizeof(double));
2613  dmtotal_Gfactor = (double*) malloc(size * sizeof(double));
2614 
2615  dmf = (double*) malloc(size * sizeof(double));
2616  acu_dmf = (double*) malloc(size * sizeof(double));
2617  inj_velocity = (double*) malloc(size * sizeof(double));
2618  virtual_velocity = (double*) malloc(size * sizeof(double));
2619 
2620  HRF = (double*) malloc(size * sizeof(double));
2621  HRF_PMX = (double*) malloc(size * sizeof(double));
2622  ROHR = (double*) malloc(size * sizeof(double));
2623  T_cyl = (double*) malloc(size * sizeof(double));
2624  acu_Mbb = (double*) malloc(size * sizeof(double));
2625  V_cyl = (double*) malloc(size * sizeof(double));
2626  Mbb = (double*) malloc(size * sizeof(double));
2627  Yair = (double*) malloc(size * sizeof(double));
2628  Yburned = (double*) malloc(size * sizeof(double));
2629  Yfuel = (double*) malloc(size * sizeof(double));
2630  U = (double*) malloc(size * sizeof(double));
2631  CV = (double*) malloc(size * sizeof(double));
2632  H_cooler = (double*) malloc(size * sizeof(double));
2633  Qcylhead = (double*) malloc(size * sizeof(double));
2634  Qcyl = (double*) malloc(size * sizeof(double));
2635  Qpis = (double*) malloc(size * sizeof(double));
2636  H = (double*) malloc(size * sizeof(double));
2637  defor = (double*) malloc(size * sizeof(double));
2638  Rmixture = (double*) malloc(size * sizeof(double));
2639  roair = (double*) malloc(size * sizeof(double));
2640  p_cyl = (double*) malloc(size * sizeof(double));
2641  dp_da_cyl = (double*) malloc(size * sizeof(double));
2642  XLO = (double*) malloc(size * sizeof(double));
2643  PEN_min = (double*) malloc(size * sizeof(double));
2644  PEN_max = (double*) malloc(size * sizeof(double));
2645  evol_mSoot = (double*) malloc(size * sizeof(double));
2646  mSootCil = (double*) malloc(size * sizeof(double));
2647  evol_Radiacion = (double*) malloc(size * sizeof(double));
2648 
2649  aux_vector = (double*) malloc(CAI * sizeof(double));
2650 
2651  // equilibrium NOx emission behaviour map
2652 
2653  YNOeq = (double**) malloc(52 * sizeof(double*));
2654  for(counter = 0; counter < 52; counter++) {
2655  YNOeq[counter] = (double*) malloc(12 * sizeof(double));
2656  }
2657 
2658  KdYNO = (double**) malloc(52 * sizeof(double*));
2659  for(counter = 0; counter < 52; counter++) {
2660  KdYNO[counter] = (double*) malloc(12 * sizeof(double));
2661  }
2662 
2663  // INJECTION RATE CONSTRUCTION
2664 
2665  for(counter = 0; counter < NIN; counter++) {
2666  SOI_IM[counter] = SOI[counter];
2667  EOI_IM[counter] = EOI[counter];
2668  SOC_IM[counter] = EVO;
2669  }
2670  vector_to_interpolate = (double*) malloc(size_inlet_inj * sizeof(double));
2671  for(counter = 0; counter < size_inlet_inj; counter++) {
2672  vector_to_interpolate[counter] = injection_rate[counter];
2673  }
2674 
2675  FUNCTION_FOR_INTERPOLATION(dmf, time_vector, CAD_injection_rate, vector_to_interpolate, size, size_inlet_inj, speed);
2676  free(vector_to_interpolate);
2677 
2678  for(counter = 0; counter < size; counter++) {
2679  if(time_vector[counter] <= SOI_IM[0] * 60. / (360. * speed)) {
2680  dmf[counter] = 0.;
2681  }
2682  if(time_vector[counter] >= EOI_IM[NIN - 1] * 60. / (360. * speed)) {
2683  dmf[counter] = 0.;
2684  }
2685  for(inj_counter = 1; inj_counter < NIN; inj_counter++) {
2686  if((time_vector[counter] >= EOI_IM[inj_counter - 1] * 60. / (360. * speed))
2687  && (time_vector[counter] <= SOI_IM[inj_counter] * 60. / (360. * speed))) {
2688  dmf[counter] = 0.;
2689  }
2690  }
2691 
2692  if(dmf[counter] < 0) {
2693  dmf[counter] = 0.;
2694  }
2695 
2696  }
2697 
2698  // injection rate area calculation and injection rate law scale: rate shape is kept but his dimension is changed to
2699  // agree with fuel mass inlet variable
2700 
2701  rate_area = 0;
2702  for(counter = 0; counter < size - 1; counter++) {
2703  rate_area = rate_area + (time_vector[counter + 1] - time_vector[counter]) * min(dmf[counter],
2704  dmf[counter + 1]) + (time_vector[counter + 1] - time_vector[counter]) * pow(
2705  pow((dmf[counter + 1] - dmf[counter]) / 2, 2), 0.5);
2706  }
2707 
2708  if(rate_area > 0) {
2709  for(counter = 0; counter < size; counter++) {
2710  dmf[counter] = dmf[counter] * mfuel * 1e-6 / rate_area;
2711  }
2712  }
2713 
2714  inj_num = NIN;
2715 
2718  CALCULUS_OF_ACCUMULATED_INJ_RATE(acu_dmf, dmf, time_vector, size);
2719 
2722  CALCULUS_OF_VIRTUAL_VELOCITY(inj_velocity, virtual_velocity, dmf, time_vector, rofuel, dc, n_holes, nozzle_d, D, size,
2723  PI, speed, EOI_IM, inj_num, SOI_IM, Piston_D, DBowl, CTM, CAD, Kswirl);
2724 
2727  STOICHIOMETRY_CONSTANTS(HC, &Kst1, &Kst2, &Kst3, &Kst4, &Kst5, &Kst6);
2728 
2731  num_i_IM = (int*) malloc(inj_num * sizeof(int));
2732 
2733  num_i = 0; /* Initialization of number of elements to 0 */
2734  num_j = 5; /* Number of sub-elements fixed to 5 */
2735 
2736  CALCULUS_OF_NUMBER_ELEMENTS(num_i_IM, time_vector, size, speed, SOI_IM, EOI_IM, inj_num);
2737 
2739  num_i = num_i_IM[0];
2740 
2741  POI_IM = (double**) malloc(inj_num * sizeof(double*));
2742  for(counter = 0; counter < inj_num; counter++) {
2743  num_i = num_i_IM[counter];
2744  POI_IM[counter] = (double*) malloc(num_i * sizeof(double));
2745  }
2746 
2747  POC_IM = (double**) malloc(inj_num * sizeof(double*));
2748  for(counter = 0; counter < inj_num; counter++) {
2749  element_value = num_i_IM[counter] * num_j;
2750  POC_IM[counter] = (double*) malloc(element_value * sizeof(double));
2751  }
2752 
2753  mfuel_i_IM = (double**) malloc(inj_num * sizeof(double*));
2754  for(counter = 0; counter < inj_num; counter++) {
2755  num_i = num_i_IM[counter];
2756  mfuel_i_IM[counter] = (double*) malloc(num_i * sizeof(double));
2757  }
2758 
2759  mfuel_ij_IM = (double**) malloc(inj_num * sizeof(double*));
2760  for(counter = 0; counter < inj_num; counter++) {
2761  num_i = num_i_IM[counter];
2762  mfuel_ij_IM[counter] = (double*) malloc(num_i * sizeof(double));
2763  }
2764 
2765  elementcontrol = (struct stControlElementComposition * *) malloc(15 * sizeof(struct stControlElementComposition));
2766  for(counter = 0; counter < 15; counter++) {
2767  elementcontrol[counter] = (struct stControlElementComposition*) malloc(size * sizeof(
2768  struct stControlElementComposition));
2769  }
2770 
2771  CALCULUS_OF_POI(POI_IM, mfuel_i_IM, mfuel_ij_IM, acu_dmf, time_vector, size, speed, num_i_IM, num_j, SOI_IM, EOI_IM,
2772  inj_num, elementcontrol);
2773 
2774  mixture_correction = (double*) malloc(num_j * sizeof(double));
2775 
2776  element = (struct stPropertiesElement * *) malloc(inj_num * sizeof(struct stPropertiesElement*));
2777  for(counter = 0; counter < inj_num; counter++) {
2778  element_value = num_i_IM[counter] * num_j;
2779  element[counter] = (struct stPropertiesElement*) malloc(element_value * sizeof(struct stPropertiesElement));
2780  }
2781 
2782  realelement = (struct stRealElementComposition * *) malloc(inj_num * sizeof(struct stRealElementComposition*));
2783  for(counter = 0; counter < inj_num; counter++) {
2784  element_value = num_i_IM[counter] * num_j;
2785  realelement[counter] = (struct stRealElementComposition*) malloc(element_value * sizeof(
2786  struct stRealElementComposition));
2787  }
2788 
2789  mixture_correction[0] = 0.319;
2790  mixture_correction[1] = 0.668;
2791  mixture_correction[2] = 1.;
2792  mixture_correction[3] = 1.397;
2793  mixture_correction[4] = 2.07;
2794 
2795  mf_burned = 0.;
2796  mf_burned_pmx = 0.;
2797 
2798  // inlet valve closing composition calculation
2799 
2800  mtotal_IVC = mairIVC;
2801  mEGR = mairIVC * rateBG;
2802 
2803  YO2IVC = YO2aire * (1 - rateBG * (mfuel * 1.e-3 + mfuel * 1e-3 * AFe) / (mairIVC * (1. - rateBG) + mfuel * 0.001));
2804  YCO2IVC = mmCO2 / mmfuel * mfuel * 0.001 * rateBG / (mfuel * 0.001 + mairIVC * (1 - rateBG));
2805  YH2OIVC = HC / 2. * mmH2O / mmfuel * mfuel * 0.001 * rateBG / (mfuel * 0.001 + mairIVC * (1 - rateBG));
2806  YN2IVC = 1 - YO2IVC - YCO2IVC - YH2OIVC;
2807 
2808  mO2_IVC = mairIVC * YO2IVC;
2809  mN2_IVC = mairIVC * YN2IVC;
2810  mCO2_IVC = mairIVC * YCO2IVC;
2811  mH2O_IVC = mairIVC * YH2OIVC;
2812 
2813  NOx_IVC = test_variables[14] * 1e-6 * 1.587 * mairIVC;
2814  mSOOT_IVC_A = test_variables[15] * 1e-6 * mairIVC;
2815  mSOOT_IVC_B = test_variables[15] * 1e-6 * mairIVC;
2816  mSOOT_IVC_C = test_variables[15] * 1e-6 * mairIVC;
2817  CO_IVC = 0;
2818  HC_IVC = 0;
2819 
2820  Cylinder_capacity = PI * Piston_D * Piston_D / 4. * S;
2821  VTDC = Cylinder_capacity / (GCRatio - 1.);
2822  f = mairIVC * 1000. / mfuel;
2823 
2824  // vectors initiate
2825 
2826  for(counter = 0; counter < size; counter++) {
2827 
2828  mtotal_bowl[counter] = mtotal_IVC;
2829  mO2_bowl[counter] = mO2_IVC;
2830  mN2_bowl[counter] = mN2_IVC;
2831  mCO2_bowl[counter] = mCO2_IVC;
2832  mH2O_bowl[counter] = mH2O_IVC;
2833  mHC_bowl[counter] = 0;
2834  Tsq_cyl[counter] = 0;
2835  Tadib_cyl[counter] = 0;
2836  mNOx_bowl[counter] = NOx_IVC;
2837  mCO_bowl[counter] = 0;
2838  HC_bowl[counter] = 0;
2839  mSOOT_bowl_A[counter] = mSOOT_IVC_A;
2840  mSOOT_bowl_B[counter] = mSOOT_IVC_B;
2841  mSOOT_bowl_C[counter] = mSOOT_IVC_C;
2842  HRF[counter] = 0;
2843  ROHR[counter] = 0;
2844  T_cyl[counter] = 0;
2845  acu_Mbb[counter] = 0;
2846  V_cyl[counter] = 0;
2847  Mbb[counter] = 0;
2848  Yair[counter] = 0;
2849  Yburned[counter] = 0;
2850  Yfuel[counter] = 0;
2851  U[counter] = 0;
2852  CV[counter] = 0;
2853  H_cooler[counter] = 0;
2854  H[counter] = 0;
2855  defor[counter] = 0;
2856  Rmixture[counter] = 0;
2857  roair[counter] = 0;
2858  p_cyl[counter] = 0;
2859  mtotal_Vc[counter] = 0;
2860  mO2_Vc[counter] = 0;
2861  mN2_Vc[counter] = 0;
2862  mCO2_Vc[counter] = 0;
2863  mH2O_Vc[counter] = 0;
2864  mNOx_Vc[counter] = 0;
2865  mSOOT_Vc_A[counter] = 0;
2866  mSOOT_Vc_B[counter] = 0;
2867  mHC_Vc[counter] = 0;
2868  dmtotal_Gfactor[counter] = 0;
2869  Qcyl[counter] = 0;
2870  Qcylhead[counter] = 0;
2871  Qpis[counter] = 0;
2872  evol_mSoot[counter] = 0.;
2873  evol_Radiacion[counter] = 0.;
2874  }
2875 
2876  p_cyl[0] = 0;
2877  p_cyl[1] = 0;
2878  T_cyl[0] = temperatureIVC;
2879  T_cyl[1] = temperatureIVC;
2880  Mbb[0] = 0;
2881  Mbb[1] = 0;
2882 
2883  // element matrix initiate: m counter indicates injection number, i counter indicates element number, j counter indicates sub-element number
2884 
2885  for(m = 0; m < inj_num; m++) {
2886  for(i = 0; i < num_i_IM[m]; i++) {
2887  for(j = 0; j < num_j; j++) {
2888 
2889  aux = i * num_j + j;
2890  element[m][aux].state = stEvaporated;
2891  element[m][aux].mtotal = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2892  element[m][aux].mf_jet = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2893  element[m][aux].mf_reac = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2894  element[m][aux].mf_evap = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2895 
2896  element[m][aux].mO2 = 0;
2897  element[m][aux].mN2 = 0;
2898  element[m][aux].mCO2 = 0.;
2899  element[m][aux].mH2O = 0.;
2900 
2901  element[m][aux].dmtotal = 0;
2902  POC_IM[m][aux] = 8888;
2903  element[m][aux].RID = 0;
2904  element[m][aux].Rpmx = 0;
2905  element[m][aux].Rpmx_value = 0;
2906  element[m][aux].mNOx = 0;
2907  element[m][aux].dNOx = 0;
2908  element[m][aux].CO = 0;
2909  element[m][aux].HC = 0;
2910  element[m][aux].mSOOT_A = 0;
2911  element[m][aux].dSOOT_A = 0;
2912  element[m][aux].mSOOT_B = 0;
2913  element[m][aux].dSOOT_B = 0;
2914  element[m][aux].mSOOT_C = 0;
2915  element[m][aux].dSOOT_C = 0;
2916  element[m][aux].TSD = 0;
2917  element[m][aux].Tadib = 0.;
2918  element[m][aux].Pcyl_POC = 1;
2919  element[m][aux].TNOx = 0;
2920  element[m][aux].X = 0;
2921  element[m][aux].soot_precursor = 0;
2922  element[m][aux].C = 1;
2923  element[m][aux].FI = 0;
2924  element[m][aux].FRLOL = 0;
2925  element[m][aux].tLOL =
2926  -11.; // This initial value serves as a way to know whether the element has reached the LOL or not.
2927 
2928  realelement[m][aux].mtotal = 0;
2929  realelement[m][aux].mO2 = 0;
2930  realelement[m][aux].mO = 0;
2931  realelement[m][aux].mCH = 0;
2932  realelement[m][aux].mN2 = 0;
2933  realelement[m][aux].mCO2 = 0;
2934  realelement[m][aux].mH2O = 0;
2935  realelement[m][aux].mf_reac = 0;
2936  realelement[m][aux].dmf_reac = 0;
2937  realelement[m][aux].mf_reac_pmx = 0;
2938 
2939  }
2940  }
2941  }
2942  for(i = 0; i < 15; i++) {
2943  for(j = 0; j < size; j++) {
2944  elementcontrol[i][j].mfuel = 0;
2945  elementcontrol[i][j].mfuel_real = 0;
2946  elementcontrol[i][j].mtotal = 0;
2947  elementcontrol[i][j].mO2 = 0;
2948  elementcontrol[i][j].mO2_real = 0;
2949  elementcontrol[i][j].mO = 0;
2950  elementcontrol[i][j].mCO2 = 0;
2951  elementcontrol[i][j].mH2O = 0;
2952  elementcontrol[i][j].mCH = 0;
2953  elementcontrol[i][j].TSD = 0;
2954  elementcontrol[i][j].Tadib = 0.;
2955  elementcontrol[i][j].dNOx = 0;
2956  elementcontrol[i][j].mNOx = 0;
2957  elementcontrol[i][j].mSOOT_A = 0;
2958  elementcontrol[i][j].mSOOT_B = 0;
2959  elementcontrol[i][j].mSOOT_C = 0;
2960  elementcontrol[i][j].X = 0;
2961  elementcontrol[i][j].FI = 0;
2962 
2963  }
2964  }
2965 
2966  // Abro el fichero interno
2967  FILE *finterno;
2968  if(RadCalc == 1 && Ang_Grab > -179.) {
2969  finterno = fopen("paquete.csv", "w");
2970  if(finterno == NULL) {
2971  printf("Error abriendo fichero interno.");
2972  // getch();
2973  exit(-1);
2974  }
2975  fprintf(finterno, "Ang,m_Soot,dm_Soot,mf,dmf \n");
2976  }
2977 
2978  FRLOL = 0.;
2979 
2982  for(counter = 0; counter < size; counter++) {
2983 
2984  mf_burned = 0;
2985  mf_burned_pmx = 0;
2986 
2987  // function call to calculate instantaneous variables with HRF
2988  CALCULATE_CYCLE(roair, CAD, delta_t, V_cyl, VTDC, counter, speed, p_cyl, HRF, acu_dmf, Mbb, acu_Mbb, AFe, f,
2989  mfuel * 1e-6, mEGR * 1e-3, mairIVC * 1e-3, T_cyl, HP, Yair, Yfuel, Yburned, U, CV,
2990  H_cooler, H, T_Evaporation_fuel, inj_fuel_temp, PRECISION_ITERATION, defor, Rmixture, Atmosphere_press, &Gamma, PI,
2991  Runiv, Piston_D, S, Crank_L, Connecting_Rod_L, E, Piston_Axis_D,
2992  Piston_Crown_H, DBowl, VBowl, M_Connecting_Rod, M_P_R_PA, MW_air, MW_fuel, MW_burned, C_ESteel, C_Mech_Defor, CTM, WC1A,
2993  WC1B, C2, C_MBLBY, Cbb, TPIS, TCYL_HEAD, TCYL, Qcylhead, Qcyl,
2994  Qpis);
2995 
2996  // mass exchange between bowl and dead volumen with dmtotal_Gfactor quantity.
2997 
2998  if(counter > 0) {
2999 
3000  dmtotal_Gfactor[counter] = (Vc_factor - 1) * (exp(-6.9 * pow(CAD[counter] / 180,
3001  2)) - exp(-6.9 * pow(CAD[counter - 1] / 180, 2))) * (mtotal_bowl[counter - 1]);
3002 
3003  mtotal_bowl[counter] = mtotal_bowl[counter - 1] + dmtotal_Gfactor[counter];
3004  mtotal_Vc[counter] = mtotal_Vc[counter - 1] - dmtotal_Gfactor[counter];
3005 
3006  if(dmtotal_Gfactor[counter] < 0) {
3007 
3008  mO2_bowl[counter] = mO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mO2_bowl[counter - 1] / mtotal_bowl[counter -
3009  1]);
3010  mN2_bowl[counter] = mN2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mN2_bowl[counter - 1] / mtotal_bowl[counter -
3011  1]);
3012  mCO2_bowl[counter] = mCO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mCO2_bowl[counter - 1] / mtotal_bowl[counter -
3013  1]);
3014  mH2O_bowl[counter] = mH2O_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mH2O_bowl[counter - 1] / mtotal_bowl[counter -
3015  1]);
3016  mNOx_bowl[counter] = mNOx_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mNOx_bowl[counter - 1] / mtotal_bowl[counter -
3017  1]);
3018  mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_bowl_A[counter - 1] /
3019  mtotal_bowl[counter - 1]);
3020  mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_bowl_B[counter - 1] /
3021  mtotal_bowl[counter - 1]);
3022  mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_bowl_C[counter - 1] /
3023  mtotal_bowl[counter - 1]);
3024 
3025  mHC_bowl[counter] = mHC_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mHC_bowl[counter - 1] / mtotal_bowl[counter -
3026  1]);
3027 
3028  mO2_Vc[counter] = mO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mO2_bowl[counter - 1] / mtotal_bowl[counter - 1]);
3029  mN2_Vc[counter] = mN2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mN2_bowl[counter - 1] / mtotal_bowl[counter - 1]);
3030  mCO2_Vc[counter] = mCO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mCO2_bowl[counter - 1] / mtotal_bowl[counter -
3031  1]);
3032  mH2O_Vc[counter] = mH2O_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mH2O_bowl[counter - 1] / mtotal_bowl[counter -
3033  1]);
3034  mNOx_Vc[counter] = mNOx_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mNOx_bowl[counter - 1] / mtotal_bowl[counter -
3035  1]);
3036  mSOOT_Vc_A[counter] = mSOOT_Vc_A[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_bowl_A[counter - 1] /
3037  mtotal_bowl[counter - 1]);
3038  mSOOT_Vc_B[counter] = mSOOT_Vc_B[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_bowl_B[counter - 1] /
3039  mtotal_bowl[counter - 1]);
3040  mSOOT_Vc_C[counter] = mSOOT_Vc_C[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_bowl_C[counter - 1] /
3041  mtotal_bowl[counter - 1]);
3042 
3043  mHC_Vc[counter] = mHC_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mHC_bowl[counter - 1] / mtotal_bowl[counter - 1]);
3044 
3045  } else {
3046 
3047  mO2_bowl[counter] = mO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mO2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3048  mN2_bowl[counter] = mN2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mN2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3049  mCO2_bowl[counter] = mCO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mCO2_Vc[counter - 1] / mtotal_Vc[counter -
3050  1]);
3051  mH2O_bowl[counter] = mH2O_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mH2O_Vc[counter - 1] / mtotal_Vc[counter -
3052  1]);
3053  mNOx_bowl[counter] = mNOx_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mNOx_Vc[counter - 1] / mtotal_Vc[counter -
3054  1]);
3055  mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_Vc_A[counter - 1] /
3056  mtotal_Vc[counter - 1]);
3057  mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_Vc_B[counter - 1] /
3058  mtotal_Vc[counter - 1]);
3059  mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_Vc_C[counter - 1] /
3060  mtotal_Vc[counter - 1]);
3061 
3062  mHC_bowl[counter] = mHC_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mHC_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3063 
3064  mO2_Vc[counter] = mO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mO2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3065  mN2_Vc[counter] = mN2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mN2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3066  mCO2_Vc[counter] = mCO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mCO2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3067  mH2O_Vc[counter] = mH2O_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mH2O_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3068  mNOx_Vc[counter] = mNOx_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mNOx_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3069  mSOOT_Vc_A[counter] = mSOOT_Vc_A[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_Vc_A[counter - 1] / mtotal_Vc[counter
3070  - 1]);
3071  mSOOT_Vc_B[counter] = mSOOT_Vc_B[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_Vc_B[counter - 1] / mtotal_Vc[counter
3072  - 1]);
3073  mSOOT_Vc_C[counter] = mSOOT_Vc_C[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_Vc_C[counter - 1] / mtotal_Vc[counter
3074  - 1]);
3075 
3076  mHC_Vc[counter] = mHC_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mHC_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3077  }
3078  } // end if counter>0
3079 
3080  // unburned in-cylinder temperrature calculation and adiabatic in-cylinder temperature calculation
3081 
3082  // TEMPERATURA SIN QUEMAR
3083  if(time_vector[counter] < (SOC_IM[0] / (6. * speed))) {
3084  Tsq_cyl[counter] = T_cyl[counter];
3085  } else {
3086  Tsq_cyl[counter] = Tsq_cyl[counter - 1] * pow((p_cyl[counter] / p_cyl[counter - 1]), ((Gamma - 1.) / Gamma));
3087  }
3088 
3089  // TEMPERATURA ADIABATICA
3090  // Tadib_cyl[counter]=Tsq_cyl[counter]+7424.7*pow(YO2IVC,0.829961)*pow((p_cyl[counter]*(1e-5)),-0.0236771);
3091 
3092  YO2i = mO2_bowl[counter] / mtotal_bowl[counter];
3093 
3094  // Version 1
3095  // Tadib_cyl[counter]=T_cyl[counter]+37630.5*YO2i/(YO2i+AFe*YO2aire);
3096 
3097  // Version 2
3098  Tadib_cyl[counter] = Tsq_cyl[counter] + 37630.5 * YO2IVC / (YO2IVC + AFe * YO2aire);
3099 
3100  if(Tadib_cyl[counter] < 2600) {
3101  Tadib_cyl[counter] = Tadib_cyl[counter] - 1.5538 * 1e-7 * pow(Tadib_cyl[counter], 2.6774);
3102  } else {
3103  Tadib_cyl[counter] = Tadib_cyl[counter] - 7.136 * 1e-10 * pow(Tadib_cyl[counter], 3.3596);
3104  }
3105 
3106  // Aminoramos esta temperatura adiabatica (por transf. de calor) segun se ha comprobado.
3107  Tadib_cyl[counter] *= 0.9;
3108 
3109  // PARAMETROS PARA CALCULO DE TEMPERATURA (RELACIONES DE ESTADO)
3110  Zst = YO2i / (YO2i + AFe * YO2aire);
3111  Z0 = -T_cyl[counter] * Zst / (Tadib_cyl[counter] - T_cyl[counter]);
3112  Z1 = 1. + inj_fuel_temp * (1. - Zst) / (Tadib_cyl[counter] - inj_fuel_temp);
3113 
3114  // lift off calculation
3115  double Rhoa = 0., Ta = 0., d0 = 0., EfectoO2 = 0., Pres = 0.,
3116  valor = 0.; // valor se usa mas abajo = 0., en el calculo de Fsoot.
3117  double x1 = 0., x2 = 0., t1 = 0., t2 = 0., t0 = 0.,
3118  K = 0.; // Parametros para calcular con exactitud las condiciones en el LOL.
3119 
3120  YO2_bowl_i = YO2IVC; // It considers YO2 at the beginning
3121  Uo_i = virtual_velocity[counter];
3122  Pres = p_cyl[counter];
3123  Ta = Tsq_cyl[counter];
3124  Rhoa = Pres / (Runiv / MW_air * Ta); // Densidad de lo que se engloba.
3125  d0 = nozzle_d;
3126  EfectoO2 = (1. + 3.33 / YO2_bowl_i) / (1. + 3.33 / YO2aire);
3127  XLO[counter] = 1.15e11 * Uo_i * pow(Ta, -3.74) * pow(Rhoa, -0.85) * pow(d0 * 1000,
3128  0.34) * EfectoO2 / 1000. * 0.62 / Kmixture1;
3129  // en metros!!
3130  PEN_min[counter] = XLO[counter];
3131  PEN_max[counter] = 0.;
3132 
3133  // mixture model begin if the injection has begun
3134 
3135  if(time_vector[counter] >= SOI_IM[0] * 60. / (360. * speed)) {
3136 
3137  // m determines the injection number; i determines the element number and j determines the sub-element number (the trajectory)
3138 
3139  for(m = 0; m < inj_num; m++) {
3140  for(i = 0; i < num_i_IM[m]; i++) {
3141  for(j = 0; j < num_j; j++) {
3142 
3143  // aux converts the element and the sub-element position into only one position
3144 
3145  aux = i * num_j + j;
3146 
3147  if((POI_IM[m][i] <= time_vector[counter]) && (element[m][aux].mf_jet > 0) && (mO2_bowl[counter] > 0)
3148  && (element[m][aux].state != stOvermixed)) {
3149 
3150  // mixture process and air mass include into element
3151 
3152  element[m][aux].dmtotal = Kmixture1 * (pow(element[m][aux].mf_jet,
3153  2) / element[m][aux].mtotal) * virtual_velocity[counter] * pow(roair[counter], 0.5) * pow(
3154  nozzle_d * 1000, -1) * pow((mO2_bowl[counter] / mtotal_bowl[counter]),
3155  Kmixture2) * (1 / mixture_correction[j]) * delta_t;
3156 
3157  element[m][aux].mtotal = element[m][aux].mtotal + element[m][aux].dmtotal;
3158  element[m][aux].mO2 = element[m][aux].mO2 + element[m][aux].dmtotal * (mO2_bowl[counter] / mtotal_bowl[counter]);
3159  element[m][aux].mCO2 = element[m][aux].mCO2 + element[m][aux].dmtotal * (mCO2_bowl[counter] / mtotal_bowl[counter]);
3160  element[m][aux].mH2O = element[m][aux].mH2O + element[m][aux].dmtotal * (mH2O_bowl[counter] / mtotal_bowl[counter]);
3161  element[m][aux].mN2 = element[m][aux].mN2 + element[m][aux].dmtotal * (mN2_bowl[counter] / mtotal_bowl[counter]);
3162  element[m][aux].mf_reac = element[m][aux].mf_reac + element[m][aux].dmtotal * (mHC_bowl[counter] /
3163  mtotal_bowl[counter]);
3164  element[m][aux].mf_evap = element[m][aux].mf_evap + element[m][aux].dmtotal * (mHC_bowl[counter] /
3165  mtotal_bowl[counter]);
3166 
3167  // element parameters to determinate combustion process
3168 
3169  element[m][aux].FI = element[m][aux].mf_reac * Kst1 / element[m][aux].mO2;
3170  element[m][aux].C = element[m][aux].mf_reac / element[m][aux].mtotal;
3171 
3172  // element penetration length
3173 
3174  element[m][aux].X += element[m][aux].C * virtual_velocity[counter] * delta_t;
3175 
3176  // Checks if x > LOL
3177  if(element[m][aux].tLOL < -10) {
3178  // In this case the element didn't reached LOL.
3179  if(element[m][aux].X >= XLO[counter]) {
3180  // Checks if it reaches LOL now.
3181  x2 = element[m][aux].X;
3182  x1 = x2 - element[m][aux].C * virtual_velocity[counter] * delta_t;
3183  if(x2 != x1) {
3184  t1 = time_vector[counter - 1];
3185  t2 = time_vector[counter];
3186  t0 = (x2 * x2 * t1 - x1 * x1 * t2) / (x2 * x2 - x1 * x1);
3187  K = x2 / sqrt(t2 - t0);
3188 
3189  element[m][aux].tLOL = t0 + pow(XLO[counter] / K, 2);
3190 
3191  // C en LOL
3192  K = element[m][aux].C * x2 / XLO[counter];
3193 
3194  // F en LOL
3195  K = 1. / (1. / K - 1.);
3196 
3197  // Fr en LOL
3198  element[m][aux].FRLOL = K * Kst1 / (mO2_bowl[counter] / mtotal_bowl[counter]);
3199  FRLOL = element[m][aux].FRLOL;
3200  } else {
3201  element[m][aux].tLOL = time_vector[counter];
3202  element[m][aux].FRLOL = element[m][aux].FI;
3203  FRLOL = element[m][aux].FRLOL;
3204  }
3205  }
3206  }
3207 
3208  // NOx mass included into element
3209 
3210  element[m][aux].mNOx = element[m][aux].mNOx + element[m][aux].dmtotal * (mNOx_bowl[counter] / mtotal_bowl[counter]);
3211 
3212  // SOOT mass of the IVC is included into element and it is oxidated totally is it is inside the flame
3213 
3214  element[m][aux].dSOOT_A = element[m][aux].dmtotal * (mSOOT_bowl_A[counter] / mtotal_bowl[counter]);
3215  element[m][aux].dSOOT_B = element[m][aux].dmtotal * (mSOOT_bowl_B[counter] / mtotal_bowl[counter]);
3216  element[m][aux].dSOOT_C = element[m][aux].dmtotal * (mSOOT_bowl_C[counter] / mtotal_bowl[counter]);
3217 
3218  // soot model A
3219  if((time_vector[counter] < (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3220  element[m][aux].mSOOT_A = element[m][aux].mSOOT_A + element[m][aux].dSOOT_A;
3221  mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter] - element[m][aux].dSOOT_A;
3222  }
3223  if((time_vector[counter] > (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3224  mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter] - element[m][aux].dSOOT_A;
3225  }
3226 
3227  // soot model B
3228 
3229  if((time_vector[counter] < (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3230  element[m][aux].mSOOT_B = element[m][aux].mSOOT_B + element[m][aux].dSOOT_B;
3231  mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter] - element[m][aux].dSOOT_B;
3232  }
3233  if((time_vector[counter] > (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3234  mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter] - element[m][aux].dSOOT_B;
3235  }
3236 
3237  // soot model C
3238 
3239  if((time_vector[counter] < (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3240  element[m][aux].mSOOT_C = element[m][aux].mSOOT_C + element[m][aux].dSOOT_C;
3241  mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter] - element[m][aux].dSOOT_C;
3242  }
3243  if((time_vector[counter] > (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3244  mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter] - element[m][aux].dSOOT_C;
3245  }
3246 
3247  // ignition delay determination
3248 
3249  if((element[m][aux].state == stEvaporated) && (SOC_IM[m] == EVO)) {
3250 
3251  Y1 = -KID3 * pow(element[m][aux].mf_reac / element[m][aux].mtotal,
3252  KID4) * pow(element[m][aux].mO2 / element[m][aux].mtotal, KID5);
3253 
3254  if(element[m][aux].mf_reac / element[m][aux].mtotal == 1) {
3255  T1 = 1000. / 353.;
3256  } else {
3257  T1 = 1000. / (T_cyl[counter] * (1 - element[m][aux].mf_reac / element[m][aux].mtotal));
3258  }
3259 
3260  if(T1 >= 2) {
3261  RID1 = Y1 + 0.155 * T1 - 0.2175;
3262  } else if(((T1 - 0.35) > 0.8) && (T1 < 2)) {
3263  RID1 = (Y1 + pow(exp((T1 - 0.35) / 5), 7) - 3.1);
3264  element[m][aux].RID = element[m][aux].RID + KID1 * pow(p_cyl[counter] * 1e-5, KID2) * pow(10, -RID1) * delta_t;
3265  } else {
3266  RID1 = (Y1 + pow(log(T1 - 0.35), 3));
3267  element[m][aux].RID = element[m][aux].RID + KID1 * pow(p_cyl[counter] * 1e-5, KID2) * pow(10, -RID1) * delta_t;
3268  }
3269  }
3270 
3271  if((element[m][aux].RID >= 1) && (SOC_IM[m] == EVO)) {
3272  SOC_IM[m] = time_vector[counter] * 6 * speed;
3273  // Se activa el calculo de la radiacion para ese pulso:
3274  CalcRad[m] = 1;
3275  }
3276 
3279  if((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3280 
3281  element[m][aux].Rpmx = element[m][aux].Rpmx + Kpmx1 * pow((mO2_bowl[counter] / mtotal_bowl[counter]),
3282  0.855) * ((1. / 4. - 1.) / pow(log(.1), 2) * pow(
3283  log(element[m][aux].FI), 2) + 1.) * exp(-Kpmx2 / T_cyl[counter]) * pow(roair[counter], Kpmx5) * pow(speed / 1000.,
3284  Kpmx3) * pow(PCR / 500., Kpmx4) * delta_t;
3285  // element[m][aux].Rpmx=element[m][aux].Rpmx+Kpmx1*pow((mO2_bowl[counter]/mtotal_bowl[counter]),Kpmx3)*((1./Kpmx4-1.)/pow(log(.1),2)*pow(log(element[m][aux].FI),2)+1.)*exp(-Kpmx2/T_cyl[counter])*pow(roair[counter],Kpmx5)*delta_t;
3286  // element[m][aux].Rpmx=element[m][aux].Rpmx+Kpmx1*pow((mO2_bowl[counter]/mtotal_bowl[counter]),Kpmx3)*pow(fabs(log(element[m][aux].FI)),Kpmx4)*exp(-Kpmx2/T_cyl[counter])*pow(roair[counter],Kpmx5)*delta_t;
3287  // element[m][aux].Rpmx_value=Kpmx1*pow((mO2_bowl[counter]/mtotal_bowl[counter]),Kpmx3)*pow(fabs(log(element[m][aux].FI)),Kpmx4)*exp(-Kpmx2/T_cyl[counter])*pow(roair[counter],Kpmx5)*delta_t;
3288  if(element[m][aux].Rpmx > 1.)
3289  element[m][aux].Rpmx = 1.;
3290  // Limitamos a 1. el valor
3291 
3292  if(element[m][aux].Pcyl_POC < 2.)
3293  // En este caso no esta inicializado: lo inicializo ahora por si acaso.
3294  element[m][aux].Pcyl_POC = p_cyl[counter];
3295 
3296  }
3297 
3298  // PREMIX COMBUSTION---------------------------------------------------------------
3299 
3300  if((element[m][aux].state == stEvaporated) && (POI_IM[m][i] < SOC_IM[m] / (6 * speed))
3301  && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3302  if(POI_IM[m][i] < SOC_IM[m] / (6 * speed)) {
3303  // Determinacion del POC (es en el inicio de la combustion premezclada)
3304  POC_IM[m][aux] = time_vector[counter];
3305 
3306  if(element[m][aux].FI > 1.) {
3307  // Caso Premezcla Rica
3308  realelement[m][aux].mf_reac_pmx = element[m][aux].mf_reac * (1. - 1. / element[m][aux].FI * element[m][aux].Rpmx);
3309  if(element[m][aux].Rpmx == 1.)
3310  // Damos el paquete por quemado en premezcla Rica:
3311  element[m][aux].state = stBurned_rich_premix;
3312  } else {
3313  // Caso Premezcla Pobre
3314  element[m][aux].Pcyl_POC = p_cyl[counter];
3315  realelement[m][aux].mf_reac_pmx = element[m][aux].mf_reac * (1. - element[m][aux].Rpmx);
3316  if(element[m][aux].Rpmx == 1.)
3317  // Damos el paquete por quemado en premezcla Pobre:
3318  element[m][aux].state = stBurned_poor_premix;
3319 
3320  }
3321  }
3322  }
3323 
3324  // FIRST DIFFUSION
3325  if((element[m][aux].state == stEvaporated) && (POI_IM[m][i] >= SOC_IM[m] / (6 * speed))) {
3326  element[m][aux].state = stInto_diffusion_flame;
3327 
3328  /* FUNCTION_SOOT_C(&soot_pre,element[m][aux].FI);
3329  element[m][aux].soot_precursor=soot_pre; */
3330  }
3331 
3332  if((element[m][aux].state == stInto_diffusion_flame) && (POI_IM[m][i] >= SOC_IM[m] / (6 * speed))
3333  && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3334 
3335  if(element[m][aux].FI <= 1) {
3336  POC_IM[m][aux] = time_vector[counter];
3337  element[m][aux].Pcyl_POC = p_cyl[counter];
3338  element[m][aux].state = stBurned_by_diffusion;
3339  }
3340 
3341  }
3342 
3343  // SECOND DIFFUSION
3344 
3345  if((element[m][aux].state == stBurned_rich_premix)) {
3346 
3347  if(element[m][aux].FI <= 1) {
3348  element[m][aux].state = stBurned_by_second_diffusion;
3349  POC_IM[m][aux] = time_vector[counter];
3350  element[m][aux].Pcyl_POC = p_cyl[counter];
3351  }
3352 
3353  }
3354 
3355  // real composition mass determination
3356 
3357  if((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3358  if(element[m][aux].FI > 1.) {
3359  // Paquete que esta a medio quemar, en condiciones ricas...
3360  realelement[m][aux].mO2 = element[m][aux].mO2 * (1. - element[m][aux].Rpmx);
3361  mf_old = realelement[m][aux].mf_reac;
3362  realelement[m][aux].mf_reac = element[m][aux].mf_reac * (1. - 1. / element[m][aux].FI * element[m][aux].Rpmx);
3363  realelement[m][aux].dmf_reac = realelement[m][aux].mf_reac - mf_old;
3364  } else {
3365  // Paquete que esta a medio quemar, en condiciones pobres...
3366  realelement[m][aux].mO2 = element[m][aux].mO2 - element[m][aux].mf_reac * Kst1 * element[m][aux].Rpmx;
3367  mf_old = realelement[m][aux].mf_reac;
3368  realelement[m][aux].mf_reac = element[m][aux].mf_reac * (1. - element[m][aux].Rpmx);
3369  realelement[m][aux].dmf_reac = realelement[m][aux].mf_reac - mf_old;
3370  }
3371  }
3372 
3373  if((element[m][aux].state == stBurned_rich_premix) || (element[m][aux].state == stInto_diffusion_flame)) {
3374  realelement[m][aux].mO2 = 0;
3375  mf_old = realelement[m][aux].mf_reac;
3376  realelement[m][aux].mf_reac = element[m][aux].mf_reac * (1. - 1. / element[m][aux].FI);
3377  realelement[m][aux].dmf_reac = realelement[m][aux].mf_reac - mf_old;
3378  }
3379 
3380  if((element[m][aux].state == stBurned_poor_premix) || (element[m][aux].state == stBurned_by_diffusion)
3381  || (element[m][aux].state == stBurned_by_second_diffusion)) {
3382  realelement[m][aux].mO2 = element[m][aux].mO2 - element[m][aux].mf_reac * Kst1;
3383  mf_old = realelement[m][aux].mf_reac;
3384  realelement[m][aux].mf_reac = 0;
3385  realelement[m][aux].dmf_reac = -mf_old;
3386  }
3387 
3388  // element temperature calculation
3389 
3390  // Wrong version element[m][aux].TSD=(1/(element[m][aux].mtotal-realelement[m][aux].mf_reac+2*realelement[m][aux].mf_reac))*((element[m][aux].mtotal-realelement[m][aux].mf_reac)*(T_cyl[counter]+37630.5*(element[m][aux].mf_reac-realelement[m][aux].mf_reac)/element[m][aux].mtotal)+inj_fuel_temp*(2*realelement[m][aux].mf_reac));
3391  // Corrected version element[m][aux].TSD=1./element[m][aux].mtotal*((element[m][aux].mtotal-realelement[m][aux].mf_reac)*(T_cyl[counter]+37630.5*(element[m][aux].mf_reac-realelement[m][aux].mf_reac)/element[m][aux].mtotal)+inj_fuel_temp*realelement[m][aux].mf_reac);
3392 
3393  // T con combustion, TCC
3394 
3395  if(element[m][aux].FI >= 1.) {
3396  // Dosado rico.
3397  TCC = Tadib_cyl[counter] * (Z1 - element[m][aux].C) / (Z1 - Zst);
3398  } else {
3399  // Dosado pobre.
3400  TCC = Tadib_cyl[counter] * (element[m][aux].C - Z0) / (Zst - Z0);
3401  }
3402 
3403  // T sin combustion, TSC
3404 
3405  TSC = (T_cyl[counter] * (1. - element[m][aux].C) + inj_fuel_temp * element[m][aux].C);
3406 
3407  // La temperatura real depende de si estamos antes o despues del SOC:
3408 
3409  if(time_vector[counter] > (SOC_IM[m] / (6 * speed)))
3410  // En este caso ya ha empezado la combustion
3411  element[m][aux].Tadib = TCC;
3412  else
3413  // En este caso NO ha empezado la combustion
3414  element[m][aux].Tadib = TSC;
3415 
3416  // Temperatura para el calculo de NOX
3417  element[m][aux].TNOx = element[m][aux].Tadib;
3418 
3419  // EMISSIONS FORMATION AND OXIDATIONS
3420 
3421  // NOX emissions
3422 
3423  if((element[m][aux].state == stBurned_by_diffusion) || (element[m][aux].state == stBurned_by_second_diffusion)
3424  || (element[m][aux].state == stBurned_poor_premix)) {
3425 
3426  if((time_vector[counter] == POC_IM[m][aux]) && (element[m][aux].state != stBurned_poor_premix)) {
3427  mNOx_bowl[counter] = mNOx_bowl[counter] - element[m][aux].mNOx * EC;
3428  element[m][aux].mNOx = element[m][aux].mNOx - element[m][aux].mNOx * EC;
3429  }
3430 
3431  FUNCTION_NOX(&YNOeq_value, &KdYNO_value, YNOeq, KdYNO, element[m][aux].TNOx, realelement[m][aux].mO2,
3432  element[m][aux].mtotal);
3433 
3434  if(YNOeq_value == 0) {
3435  YNOeq_value = 1e-6;
3436  }
3437 
3438  if(YNOeq_value - element[m][aux].mNOx / element[m][aux].mtotal > 0) {
3439  element[m][aux].dNOx = (KNOx2) * element[m][aux].mtotal * KdYNO_value * (YNOeq_value - element[m][aux].mNOx /
3440  element[m][aux].mtotal) / YNOeq_value * delta_t;
3441  } else {
3442  element[m][aux].dNOx = 0;
3443  }
3444 
3445  element[m][aux].mNOx = element[m][aux].mNOx + element[m][aux].dNOx;
3446  mNOx_bowl[counter] = mNOx_bowl[counter] + element[m][aux].dNOx;
3447 
3448  }
3449 
3450  // SOOT formation and soot oxidation
3451  // SOOT model A: HIROYASU
3452  // Hiroyasu formation process
3453 
3454  element[m][aux].dSOOT_A = KSOOTA1 * realelement[m][aux].mf_reac * pow(p_cyl[counter],
3455  KSOOTA2) * exp(-KSOOTA3 / element[m][aux].Tadib) * delta_t;
3456  element[m][aux].mSOOT_A = element[m][aux].mSOOT_A + element[m][aux].dSOOT_A;
3457 
3458  // Hiroyasu oxidation process
3459 
3460  element[m][aux].dSOOT_A = KSOOTA4 * element[m][aux].mSOOT_A * pow((realelement[m][aux].mO2 / element[m][aux].mtotal),
3461  KSOOTA5) * pow(p_cyl[counter], KSOOTA6) * exp(
3462  -KSOOTA7 / element[m][aux].Tadib) * delta_t;
3463  if(element[m][aux].dSOOT_A <= element[m][aux].mSOOT_A) {
3464  element[m][aux].mSOOT_A = element[m][aux].mSOOT_A - element[m][aux].dSOOT_A;
3465  } else {
3466  element[m][aux].mSOOT_A = 0;
3467  }
3468 
3469  // SOOT model C
3470 
3471  // CMT formation-oxidation process (both processes "integrated")
3472 
3473  // Oxidacion (debido a la desaparicion de masa de combustible):
3474  if(realelement[m][aux].dmf_reac < 0.)
3475  element[m][aux].mSOOT_C *= realelement[m][aux].mf_reac / (realelement[m][aux].mf_reac - realelement[m][aux].dmf_reac);
3476 
3477  // Formacion:
3478  if(element[m][aux].tLOL > -10.) {
3479  // Checks if the element has reached the LOL
3480  // element[m][aux].mSOOT_C=KSOOTC1*realelement[m][aux].mf_reac*pow(Pres/(14.8*Runiv/MW_air*1000.),2.2)*(-0.80545455/element[m][aux].FRLOL+0.49715984)*pow((time_vector[counter]-element[m][aux].tLOL)*1.e6,0.6)*exp(-2399./Tadib_cyl[counter]);
3481  if((time_vector[counter] - element[m][aux].tLOL) > 0.) {
3482  // valor=KSOOTC1*realelement[m][aux].mf_reac*pow(Pres/(14.8*Runiv/MW_air*1000.),2.2)*0.15*(element[m][aux].FRLOL-2.)*1.e6*exp(-2399./Tadib_cyl[counter])*delta_t;
3483  valor = KSOOTC1 * realelement[m][aux].mf_reac * pow(Pres / (14.8 * Runiv / MW_air * 1000.),
3484  2.2) * (-1. / element[m][aux].FRLOL + 0.5) * 1.e6 * exp(
3485  -2399. / Tadib_cyl[counter]) * delta_t;
3486  // KSOOTC1 vale 0.9e-3 para el caso de Exp 1
3487  // He anadido un factor de correccion de la presion (instantanea).
3488  /* valor=time_vector[counter]-element[m][aux].tLOL-delta_t; // Valor aqui es un calculo intermedio...
3489  if(valor<0.)
3490  valor=0.; // KSOOTC1 vale 1.e-2 en el caso de Exp 0.6 para el tiempo de residencia
3491  valor=KSOOTC1*realelement[m][aux].mf_reac*pow(Pres/(14.8*Runiv/MW_air*1000.),2.2)*(-1./element[m][aux].FRLOL+0.5)*(pow((time_vector[counter]-element[m][aux].tLOL)*1.e6,0.6)-pow(valor*1.e6,0.6))*exp(-2399./Tadib_cyl[counter]); */
3492  } else
3493  valor = 0.;
3494  if(valor < 0.) // Limitamos a cero.
3495  valor = 0.;
3496 
3497  element[m][aux].mSOOT_C += valor;
3498 
3499  // Limitamos el valor de mSOOT al valor de mf_reac
3500  if(element[m][aux].mSOOT_C > realelement[m][aux].mf_reac)
3501  element[m][aux].mSOOT_C = realelement[m][aux].mf_reac;
3502  }
3503  if(element[m][aux].mSOOT_C < 0.)
3504  // Limitamos a cero.
3505  element[m][aux].mSOOT_C = 0.;
3506 
3507  // Anadimos la masa al contador de masa de Soot:
3508  evol_mSoot[counter] += element[m][aux].mSOOT_C;
3509  // La radiacion se calculara mas adelante...
3510 
3511  // HC formation by overmixed element
3512 
3513  if((SOC_IM[m] * 60. / (360. * speed) < time_vector[counter]) && (POI_IM[m][i] < SOC_IM[m] * 60. / (360. * speed))
3514  && (element[m][aux].state == stEvaporated)) {
3515  if(element[m][aux].FI < 0.5) {
3516 
3517  mHC_bowl[counter] = mHC_bowl[counter] + element[m][aux].mf_reac;
3518  element[m][aux].state = stOvermixed;
3519  realelement[m][aux].mf_reac = 0;
3520  element[m][aux].mf_reac = 0;
3521  element[m][aux].mf_jet = 0;
3522 
3523  }
3524  }
3525  } // if time>POI
3526  } // for j
3527 
3528  // element control assignation
3529  for(i_aux = 0; i_aux < 5; i_aux++) {
3530  j_aux = i_aux * 3;
3531 
3532  if((elementcontrol[j_aux][0].inj_number == m) && (elementcontrol[j_aux][0].num_i == i)) {
3533  aux = i * num_j + 0;
3534  j_aux = i_aux * 3;
3535  elementcontrol[j_aux][counter].mtotal = element[m][aux].mtotal;
3536  elementcontrol[j_aux][counter].mfuel = element[m][aux].mf_reac;
3537  elementcontrol[j_aux][counter].mfuel_real = realelement[m][aux].mf_reac;
3538  elementcontrol[j_aux][counter].mO2 = element[m][aux].mO2;
3539  elementcontrol[j_aux][counter].mO2_real = realelement[m][aux].mO2;
3540  elementcontrol[j_aux][counter].TSD = element[m][aux].TSD;
3541  elementcontrol[j_aux][counter].Tadib = element[m][aux].Tadib;
3542  elementcontrol[j_aux][counter].X = element[m][aux].X;
3543  elementcontrol[j_aux][counter].FI = element[m][aux].FI;
3544 
3545  aux = i * num_j + 2;
3546  j_aux = i_aux * 3 + 1;
3547  elementcontrol[j_aux][counter].mtotal = element[m][aux].mtotal;
3548  elementcontrol[j_aux][counter].mfuel = element[m][aux].mf_reac;
3549  elementcontrol[j_aux][counter].mfuel_real = realelement[m][aux].mf_reac;
3550  elementcontrol[j_aux][counter].mO2 = element[m][aux].mO2;
3551  elementcontrol[j_aux][counter].mO2_real = realelement[m][aux].mO2;
3552  elementcontrol[j_aux][counter].TSD = element[m][aux].TSD;
3553  elementcontrol[j_aux][counter].Tadib = element[m][aux].Tadib;
3554  elementcontrol[j_aux][counter].X = element[m][aux].X;
3555  elementcontrol[j_aux][counter].FI = element[m][aux].FI;
3556 
3557  aux = i * num_j + 4;
3558  j_aux = i_aux * 3 + 2;
3559  elementcontrol[j_aux][counter].mtotal = element[m][aux].mtotal;
3560  elementcontrol[j_aux][counter].mfuel = element[m][aux].mf_reac;
3561  elementcontrol[j_aux][counter].mfuel_real = realelement[m][aux].mf_reac;
3562  elementcontrol[j_aux][counter].mO2 = element[m][aux].mO2;
3563  elementcontrol[j_aux][counter].mO2_real = realelement[m][aux].mO2;
3564  elementcontrol[j_aux][counter].TSD = element[m][aux].TSD;
3565  elementcontrol[j_aux][counter].Tadib = element[m][aux].Tadib;
3566  elementcontrol[j_aux][counter].X = element[m][aux].X;
3567  elementcontrol[j_aux][counter].FI = element[m][aux].FI;
3568 
3569  }
3570  }
3571 
3572  } // for i
3573  } // for m
3574 
3575  for(m = 0; m < inj_num; m++) {
3576  for(i = 0; i < num_i_IM[m]; i++) {
3577  for(j = 0; j < num_j; j++) {
3578  if(POI_IM[m][i] <= time_vector[counter]) {
3579  aux = i * num_j + j;
3580 
3581  if(((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed))))
3582  || (element[m][aux].state == stBurned_rich_premix) || (element[m][aux].state == stBurned_poor_premix)
3583  || (element[m][aux].state == stBurned_by_diffusion) || (element[m][aux].state == stBurned_by_second_diffusion)
3584  || (element[m][aux].state == stInto_diffusion_flame)) {
3585  mf_burned = mf_burned + (element[m][aux].mf_reac - realelement[m][aux].mf_reac);
3586  }
3587  // Para HRL_PMX
3588  if(((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed))))
3589  || (element[m][aux].state == stBurned_rich_premix) || (element[m][aux].state == stBurned_poor_premix)
3590  || (element[m][aux].state == stBurned_by_second_diffusion)) {
3591  mf_burned_pmx = mf_burned_pmx + (element[m][aux].mf_reac - realelement[m][aux].mf_reac_pmx);
3592  }
3593 
3594  }
3595  }
3596  }
3597  }
3598 
3599  } // if time>SOI
3600  HRF[counter] = mf_burned / (mfuel * 1.e-3);
3601  HRF_PMX[counter] = mf_burned_pmx / (mfuel * 1.e-3);
3602 
3603  if(counter > 0) {
3604  mO2_bowl[counter] = mO2_bowl[counter] - (HRF[counter] - HRF[counter - 1]) * mfuel * 1.e-3 * Kst1;
3605  mCO2_bowl[counter] = mCO2_bowl[counter] + (HRF[counter] - HRF[counter - 1]) * mfuel * 1.e-3 * Kst5;
3606  mH2O_bowl[counter] = mH2O_bowl[counter] + (HRF[counter] - HRF[counter - 1]) * mfuel * 1.e-3 * Kst6;
3607  }
3608 
3609  mSootCil[counter] = mSOOT_IVC_C * (mO2_bowl[counter] + mO2_Vc[counter]) / mO2_IVC;
3610 
3615  double lambda = 0., dlambda = 0., tgAng = 0., Teje = 0., Tllama = 0., Xeje = 0., radio = 0., aux1 = 0., Xmax = 0.;
3616  int k = 0, iini = 0, ifin = 0;
3617 
3618  if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3619  && (RadCalc == 1)) {
3620 
3621  // Cabecera del fichero 1:
3622  fprintf(foculto, "Time,%le,[s],Ang,%le,[cig],LOL,%le,[mm],FrLOL,%le\n#iny,Elem,r,x,msoot,mtot,rho,T,FI,mf_reac\n",
3623  time_vector[counter], Ang_Grab, XLO[counter] * 1000., FRLOL);
3624 
3625  // Cuerpo del fichero 1:
3626  for(m = 0; m < inj_num; m++) {
3627  for(i = 0; i < num_i_IM[m]; i++) {
3628  for(j = 0; j < num_j; j++) {
3629 
3630  // aux converts the element and the sub-element position into only one position
3631 
3632  aux = i * num_j + j;
3633  fprintf(foculto, "%d,%d,%d,%le,%le,%le,%le,%le,%le,%le\n", m + 1, i + 1, j + 1, element[m][aux].X,
3634  element[m][aux].mSOOT_C, element[m][aux].mtotal, roair[counter],
3635  element[m][aux].Tadib, element[m][aux].FI, realelement[m][aux].mf_reac);
3636  }
3637  }
3638  }
3639  }
3640 
3641  // Dimensionamos el array segun el numero de pulsos
3642  stRadArray **Radiation;
3643  Radiation = (struct stRadArray * *) malloc(inj_num * sizeof(struct stRadArray*));
3644  Itot = (double*) malloc(inj_num * sizeof(double));
3645 
3646  // Calculamos las cosas comunes. Despues ya se hara el barrido de lambda.
3647 
3648  // Lo primero es saber cuantos puntos hay (puntos del paquete 5 con Xsoot distinto de 0):
3649  for(m = 0; m < inj_num; m++) {
3650  iini = 0;
3651  ifin = 0;
3652  Itot[m] = 0; // Inicializamos a 0 la intensidad.
3653  Xmax = 0.;
3654  if(CalcRad[m] == 1 && RadCalc == 1) { // Solo se calcula si toca y si el usuario quiere
3655  j = 4; // Paquete 5
3656  i = 0;
3657  aux = i * num_j + j;
3658  // Busco el primer punto en el eje con masa de Soot distinta de 0
3659  while((element[m][aux].mSOOT_C < 1.e-20) && (i < (num_i_IM[m] - 1))) {
3660  i++;
3661  aux = i * num_j + j;
3662  }
3663  iini = i;
3664  while((element[m][aux].mSOOT_C > 0.) && (i < (num_i_IM[m] - 1))) {
3665  i++;
3666  aux = i * num_j + j;
3667  }
3668  if(i == (num_i_IM[m] - 1)) {
3669  if(element[m][aux].mSOOT_C > 0.)
3670  ifin = i;
3671  else
3672  ifin = i - 1;
3673  } else
3674  ifin = i - 1;
3675 
3676  if(iini < ifin) {
3677  // Solo en este caso conviene proseguir con el calculo...
3678  // Dimensionamos el vector. Tengo en cuenta que anadira un punto por delante (el que esta justo en el frente) y otro por detras (el que esta justo en el LOL)
3679  Radiation[m] = (struct stRadArray*) malloc((ifin - iini + 3) * sizeof(struct stRadArray));
3680 
3681  // Inicializaciones
3682  tgAng = tan(ANG_CHORRO / 2 * PI / 180.) * Kmixture1 / 0.75;
3683 
3684  // El primer elemento es el del frente de llama. Busco con precision su posicion.
3685  if(iini > 0) {
3686  // En este caso el eje del chorro ya ha llegado a estequiometrico. El punto esta entre iini e iini-1.
3687  Radiation[m][0].x = (element[m][iini * num_j + 4].X * (element[m][(iini - 1) * num_j + 4].FI - 1.) + element[m][(iini -
3688  1) * num_j + 4].X * (1. - element[m][iini * num_j + 4].FI)) / (element[m][(iini - 1) * num_j + 4].FI - element[m][iini *
3689  num_j + 4].FI);
3690  if(PEN_max[counter] < Radiation[m][0].x)
3691  PEN_max[counter] = Radiation[m][0].x;
3692  } else {
3693  // En este caso tengo que predecir la posicion de x del frente si se hubiera llegado a esa posicion:
3694  Radiation[m][0].x = element[m][4].X * Kst1 * element[m][4].mf_jet / element[m][4].mO2;
3695  if(PEN_max[counter] < element[m][4].X)
3696  PEN_max[counter] = element[m][4].X;
3697  }
3698 
3699  // Reajusto iini para asegurar utilizar el primer paquete disponible.
3700  iini--;
3701  Radiation[m][0].RFlame = 0.;
3702  Radiation[m][0].R = Radiation[m][0].x * tgAng;
3703  Radiation[m][0].np = 1;
3704  Radiation[m][0].dr = 1.e-3;
3705  Radiation[m][0].T[0] = Tadib_cyl[counter];
3706  Radiation[m][0].Xsoot[0] = 0.;
3707  Radiation[m][0].Tau[0] = 1.;
3708  Radiation[m][0].PTau[0] = 1.;
3709  Radiation[m][0].Ilambda = 0.;
3710  Radiation[m][0].Itot = 0.;
3711 
3712  for(i = (iini + 1); i <= ifin; i++) {
3713  j = 4; // Paquete 5
3714  aux = i * num_j + j;
3715  Radiation[m][i - iini].x = element[m][aux].X;
3716  aux1 = Radiation[m][i - iini].x / Radiation[m][0].x;
3717  if(aux1 < cero)
3718  aux1 = cero;
3719  aux1 = -log(aux1) / 4.6;
3720  if(aux1 < 0.)
3721  aux1 = 0.;
3722  Radiation[m][i - iini].RFlame = Radiation[m][i - iini].x * tgAng * sqrt(aux1);
3723  Radiation[m][i - iini].R = Radiation[m][i - iini].x * tgAng;
3724 
3725  if(Radiation[m][i - iini].RFlame > cero) {
3726  Radiation[m][i - iini].np = NR;
3727  // Maxima resolucion (segun tamano de vectores)
3728  Radiation[m][i - iini].dr = Radiation[m][i - iini].RFlame * 2. / (float) Radiation[m][i - iini].np;
3729  // Ahora rellenamos los vectores (en funcion de r):
3730  Teje = element[m][aux].Tadib;
3731  Tllama = Tadib_cyl[counter];
3732  Xeje = element[m][aux].mSOOT_C / element[m][aux].mtotal * p_cyl[counter] / (287. * Teje) /
3733  1800.; // Se asume que Rho_soot es 1800 kg/m^3
3734  if(Xeje > Xmax)
3735  Xmax = Xeje;
3736  for(k = 0; k < Radiation[m][i - iini].np; k++) {
3737  radio = -Radiation[m][i - iini].RFlame + (2. * k + 1.) / 2. * Radiation[m][i - iini].dr;
3738  Radiation[m][i - iini].T[k] = Teje + (Radiation[m][0].x / Radiation[m][i - iini].x * (exp(-4.6 * pow(
3739  radio / Radiation[m][i - iini].R, 2)) - 1.) / (1. - Radiation[m][0].x / Radiation[m][i - iini].x)) * (Tllama - Teje);
3740  // Radiation[m][i-iini].Xsoot[k]=Xeje*(1.-fabs(radio/Radiation[m][i-iini].RFlame));
3741  Radiation[m][i - iini].Xsoot[k] = Xeje * exp(-4.6 * pow(radio / Radiation[m][i - iini].RFlame, 2));
3742  }
3743  } else {
3744  Radiation[m][i - iini].np = 1.;
3745  Radiation[m][i - iini].dr = 0.001;
3746  Teje = element[m][aux].Tadib;
3747  if(Teje > cero)
3748  Xeje = element[m][aux].mSOOT_C / element[m][aux].mtotal * p_cyl[counter] / (287. * Teje) / 1800.;
3749  // Se asume que Rho_soot es 1800 kg/m^3
3750  if(Xeje > Xmax)
3751  Xmax = Xeje;
3752  Radiation[m][i - iini].T[0] = Teje;
3753  Radiation[m][i - iini].Xsoot[0] = Xeje;
3754  }
3755  Radiation[m][i - iini].Itot = 0.; // Se inicializa a cero la integral de la radiacion monocrometica.
3756 
3757  } // End for i
3758 
3759  // Colocamos ahora el paquete que se encuentra justo en el lift-off
3760  // Compruebo si no es el ultimo paquete
3761  if(ifin < (num_i_IM[m] - 1)) {
3762  // No es el ultimo paquete. Entonces coloco uno justo en el lift-off si no hay incoherencias...
3763  if(XLO[counter] < Radiation[m][ifin - iini].x) {
3764  Radiation[m][ifin + 1 - iini].x = XLO[counter];
3765  aux = (ifin + 1) * num_j + 4;
3766  // T se interpola
3767  Radiation[m][ifin + 1 - iini].T[0] = (element[m][aux].Tadib * (element[m][aux - num_j].X - Radiation[m][ifin + 1 -
3768  iini].x) + element[m][aux - num_j].Tadib * (Radiation[m][ifin + 1 - iini].x - element[m][aux].X)) /
3769  (element[m][aux - num_j].X - element[m][aux].X);
3770  } else {
3771  // Dado que el XLO es mayor, anado un paquete en la misma posicion
3772  Radiation[m][ifin + 1 - iini].x = Radiation[m][ifin - iini].x;
3773  Radiation[m][ifin + 1 - iini].T[0] = Radiation[m][ifin - iini].T[(Radiation[m][ifin - iini].np - 1) / 2];
3774  }
3775  } else {
3776  // Es el ultimo paquete. Coloco el paquete "extra" en el mismo sitio que el ultimo
3777  Radiation[m][ifin + 1 - iini].x = Radiation[m][ifin - iini].x;
3778  PEN_min[counter] = Radiation[m][ifin - iini].x;
3779  Radiation[m][ifin + 1 - iini].T[0] = Radiation[m][ifin - iini].T[(Radiation[m][ifin - iini].np - 1) / 2];
3780  }
3781  aux1 = Radiation[m][ifin + 1 - iini].x / Radiation[m][0].x;
3782  if(aux1 < cero)
3783  aux1 = cero;
3784  aux1 = -log(aux1) / 4.6;
3785  if(aux1 < 0.)
3786  aux1 = 0.;
3787  Radiation[m][ifin + 1 - iini].RFlame = Radiation[m][ifin + 1 - iini].x * tgAng * sqrt(aux1);
3788  Radiation[m][ifin + 1 - iini].R = Radiation[m][ifin + 1 - iini].x * tgAng;
3789  Radiation[m][ifin + 1 - iini].np = 1;
3790  Radiation[m][ifin + 1 - iini].dr = 0.001;
3791  Radiation[m][ifin + 1 - iini].Xsoot[0] = 0.;
3792  Radiation[m][ifin + 1 - iini].Tau[0] = 1.;
3793  Radiation[m][ifin + 1 - iini].PTau[0] = 1.;
3794  Radiation[m][ifin + 1 - iini].Ilambda = 0.;
3795  Radiation[m][ifin + 1 - iini].Itot = 0.;
3796 
3797  if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3798  && (RadCalc == 1)) {
3799  // Rellenamos la cabecera del fichero oculto 2:
3800  fprintf(foculto2, "Time,%le,[s],Ang,%le,[cig],LOL,%le,[mm]\n", time_vector[counter], Ang_Grab, XLO[counter] * 1000.);
3801  // Coordenada x
3802  fprintf(foculto2, "x");
3803  for(i = iini; i <= (ifin + 1); i++)
3804  fprintf(foculto2, ",%le", Radiation[m][i - iini].x);
3805  fprintf(foculto2, "\n");
3806  // Coordenada Rllama
3807  fprintf(foculto2, "Rllama");
3808  for(i = iini; i <= (ifin + 1); i++)
3809  fprintf(foculto2, ",%le", Radiation[m][i - iini].RFlame);
3810  fprintf(foculto2, "\n");
3811  // Coordenada R
3812  fprintf(foculto2, "R");
3813  for(i = iini; i <= (ifin + 1); i++)
3814  fprintf(foculto2, ",%le", Radiation[m][i - iini].R);
3815  fprintf(foculto2, "\n");
3816  // Xsoot en el eje
3817  fprintf(foculto2, "Xsoot");
3818  fprintf(foculto2, ",0.0");
3819  for(i = iini + 1; i <= ifin; i++)
3820  fprintf(foculto2, ",%le", Radiation[m][i - iini].Xsoot[(NR - 1) / 2]);
3821  fprintf(foculto2, ",0.0\n");
3822  // T en el eje
3823  fprintf(foculto2, "T");
3824  for(i = iini; i <= (ifin + 1); i++) {
3825  if(Radiation[m][i - iini].np == NR)
3826  aux1 = Radiation[m][i - iini].T[(NR - 1) / 2];
3827  else
3828  aux1 = Radiation[m][i - iini].T[0];
3829  fprintf(foculto2, ",%le", aux1);
3830  }
3831  fprintf(foculto2, "\n");
3832  }
3833 
3834  // Ahora hacemos un barrido de lambda...
3835  dlambda = 0.5;
3836  for(lambda = dlambda; lambda < 5.; lambda = lambda + dlambda) {
3837  for(i = (iini + 1); i <= ifin; i++) {
3838  // Rellenamos los vectores (en funcion de r):
3839  for(k = 0; k < Radiation[m][i - iini].np; k++) {
3840  Radiation[m][i - iini].Tau[k] = exp(-6.3 * Radiation[m][i - iini].Xsoot[k] * 1.e6 * Radiation[m][i - iini].dr / pow(
3841  lambda, 1.22 - 0.245 * log(lambda)));
3842  if(k == 0)
3843  Radiation[m][i - iini].PTau[k] = Radiation[m][i - iini].Tau[k];
3844  else
3845  Radiation[m][i - iini].PTau[k] = Radiation[m][i - iini].Tau[k] * Radiation[m][i - iini].PTau[k - 1];
3846  // Acumulo la atenuacion
3847  }
3848 
3849  // Calculamos la radiacion monocrometica:
3850  Radiation[m][i - iini].Ilambda = 0.;
3851  for(k = 0; k < Radiation[m][i - iini].np; k++) {
3852  if(Radiation[m][i - iini].T[k] > cero) {
3853  // Es el unico caso en el que tiene sentido calcularla. Si no, la radiacion es nula
3854  aux1 = 3.743e8 / (PI * pow(lambda, 5) * (exp(1.4387e4 / (lambda * Radiation[m][i - iini].T[k])) - 1.));
3855  aux1 = aux1 * (1. - Radiation[m][i - iini].Tau[k]) * Radiation[m][i - iini].PTau[k];
3856  Radiation[m][i - iini].Ilambda += aux1 * 1.e6;
3857  // Expresada en W/stdeg m^3
3858  }
3859  }
3860  Radiation[m][i - iini].Itot += Radiation[m][i - iini].Ilambda;
3861  } // End for i
3862 
3863  if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3864  && (RadCalc == 1)) {
3865  // Escribimos el resultado
3866  // Ilambda
3867  fprintf(foculto2, "%le", lambda);
3868  for(i = iini; i <= (ifin + 1); i++)
3869  fprintf(foculto2, ",%le", Radiation[m][i - iini].Ilambda);
3870  fprintf(foculto2, "\n");
3871 
3872  // KL
3873  fprintf(foculto2, "%le", lambda);
3874  fprintf(foculto2, ",0.0");
3875  for(i = iini + 1; i <= ifin; i++)
3876  fprintf(foculto2, ",%le", -log(Radiation[m][i - iini].PTau[NR - 1]));
3877  fprintf(foculto2, ",0.0\n");
3878  }
3879 
3880  } // End for lambda
3881 
3882  // Calculamos la radiacion total en cada x:
3883  for(i = (iini + 1); i <= ifin; i++) {
3884  // Descontamos la mitad del ultimo punto
3885  Radiation[m][i - iini].Itot -= Radiation[m][i - iini].Ilambda / 2.;
3886  Radiation[m][i - iini].Itot *= dlambda * 1.e-6;
3887  // Multiplicamos por dlambda.
3888  // Anadimos el tramo final (suponemos que en 10 um ya vale 0).
3889  Radiation[m][i - iini].Itot += Radiation[m][i - iini].Ilambda / 2. * (10. - lambda) * 1.e-6;
3890  // Intensidad en W/stdeg m^2
3891  }
3892 
3893  if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3894  && (RadCalc == 1)) {
3895  // Grabamos el resultado:
3896  fprintf(foculto2, "Itot");
3897  // Itot
3898  for(i = iini; i <= (ifin + 1); i++)
3899  fprintf(foculto2, ",%le", Radiation[m][i - iini].Itot);
3900  fprintf(foculto2, "\n");
3901  }
3902  // Finalmente calculamos la radiacion total:
3903  if(iini < 0) {
3904  // Se ha de quitar el primer paquete (frente), pues se trata de un chorro no estabilizado.
3905  for(i = (iini + 2); i <= (ifin + 1); i++) {
3906  Itot[m] += (Radiation[m][i - iini - 1].Itot + Radiation[m][i - iini].Itot) / 2. * (Radiation[m][i - iini - 1].x -
3907  Radiation[m][i - iini].x) * (Radiation[m][i - iini - 1].RFlame + Radiation[m][i - iini].RFlame) * PI;
3908  }
3909  } else {
3910  for(i = (iini + 1); i <= (ifin + 1); i++) {
3911  Itot[m] += (Radiation[m][i - iini - 1].Itot + Radiation[m][i - iini].Itot) / 2. * (Radiation[m][i - iini - 1].x -
3912  Radiation[m][i - iini].x) * (Radiation[m][i - iini - 1].RFlame + Radiation[m][i - iini].RFlame) * PI;
3913  }
3914  }
3915  Itot[m] *= PI;
3916  // Corregimos el valor segun la correlacion emperica obtenida.
3917  // Tengo en cuenta dos acortamientos de la llama: por YO2 < 0.23 y por chorro no estacionario.
3918  aux1 = 1. / (1. - pow(pow(Radiation[m][1].x / Radiation[m][0].x,
3919  2) * mO2_bowl[counter] / mtotal_bowl[counter] * Xmax / (0.23 * 0.004298537), 1. / 2.5));
3920  // Lo limitamos entre 1 y 2:
3921  if(aux1 > 2.)
3922  aux1 = 2.;
3923  if(aux1 < 1.)
3924  aux1 = 1.;
3925  Itot[m] *= aux1;
3926 
3927  if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3928  && (RadCalc == 1))
3929  fprintf(foculto2, "ITOT,%le\n", Itot[m]);
3930 
3931  free(Radiation[m]);
3932  } // End if iini<ifin
3933  } // End if CalcRad
3934  } // End for m
3935 
3936  // Almacenamos el valor de la intensidad de todos los pulsos en el vector de salida:
3937  evol_Radiacion[counter] = 0.;
3938  for(m = 0; m < inj_num; m++)
3939  evol_Radiacion[counter] += Itot[m];
3940  // Internamente ACT trabaja con un chorro que reune a todos los chorros (# orificios). Por tanto no hay que multiplicar por n_holes
3941  // evol_Radiacion[counter]*=n_holes;
3942 
3947  // Liberamos la memoria:
3948  free(Radiation);
3949  free(Itot);
3950 
3955  // Grabamos fichero interno
3956  elem_i = 24;
3957  pulso_i = 0;
3958  if(RadCalc == 1 && Ang_Grab > -179.)
3959  fprintf(finterno, "%le,%le,%le,%le,%le\n", time_vector[counter] * 6. * speed, element[pulso_i][elem_i].mSOOT_C,
3960  element[pulso_i][elem_i].dSOOT_C, realelement[pulso_i][elem_i].mf_reac,
3961  realelement[pulso_i][elem_i].dmf_reac);
3962 
3963  // Reajuste, por si acaso:
3964  if(PEN_min[counter] > PEN_max[counter])
3965  PEN_min[counter] = 0.;
3966 
3967  } // for counter
3968 
3969  if(RadCalc == 1 && Ang_Grab > -179.) {
3970  fclose(foculto);
3971  fclose(foculto2);
3972  }
3973 
3976  // Cierro fichero interno
3977  if(RadCalc == 1 && Ang_Grab > -179.)
3978  fclose(finterno);
3979 
3980  mSOOT_bowl_A_i_burned = mSOOT_bowl_A[size - 1];
3981  mSOOT_bowl_B_i_burned = mSOOT_bowl_B[size - 1];
3982  mSOOT_bowl_C_i_burned = mSOOT_bowl_C[size - 1];
3983 
3984  for(m = 0; m < inj_num; m++) {
3985  for(i = 0; i < num_i_IM[m]; i++) {
3986  for(j = 0; j < num_j; j++) {
3987  aux = i * num_j + j;
3988  mSOOT_bowl_A[size - 1] = mSOOT_bowl_A[size - 1] + element[m][aux].mSOOT_A;
3989  mSOOT_bowl_B[size - 1] = mSOOT_bowl_B[size - 1] + element[m][aux].mSOOT_B;
3990  mSOOT_bowl_C[size - 1] = mSOOT_bowl_C[size - 1] + element[m][aux].mSOOT_C;
3991 
3992  if(element[m][aux].FI <= 1) {
3993  mSOOT_bowl_A_i_burned = mSOOT_bowl_A_i_burned + element[m][aux].mSOOT_A;
3994  mSOOT_bowl_B_i_burned = mSOOT_bowl_B_i_burned + element[m][aux].mSOOT_B;
3995  mSOOT_bowl_C_i_burned = mSOOT_bowl_C_i_burned + element[m][aux].mSOOT_C;
3996  }
3997  }
3998  }
3999  }
4000 
4001  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4002  vector_to_interpolate[0] = 0;
4003  for(counter = 1; counter < size; counter++) {
4004  vector_to_interpolate[counter] = (HRF[counter] - HRF[counter - 1]) * (mfuel * 1e-6) * HP /
4005  (time_vector[counter] - time_vector[counter - 1]);
4006  }
4007 
4008  /* ROHR[0]=0;
4009  ROHR[size-1]=vector_to_interpolate[size-1];
4010  for(counter=1;counter<size-1;counter++){
4011  ROHR[counter]=(vector_to_interpolate[counter-1]+vector_to_interpolate[counter]+vector_to_interpolate[counter+1])/3;
4012  } */ // Quito el filtrado de la ROHR, pues ya sale bien.
4013  for(counter = 1; counter < size - 1; counter++) {
4014  ROHR[counter] = vector_to_interpolate[counter];
4015  }
4016  free(vector_to_interpolate);
4017 
4018  complete_size = 0;
4019  complete_prev_size = 0;
4020  complete_post_size = 0;
4021  counter_CAD_1 = IVC;
4022  counter_CAD_2 = EVO;
4023 
4024  do {
4025  complete_prev_size = complete_prev_size + 1;
4026  counter_CAD_1 = counter_CAD_1 - delta_t * 360. * speed / 60.;
4027  } while(counter_CAD_1 >= -180);
4028  complete_prev_size--;
4029  do {
4030  complete_post_size = complete_post_size + 1;
4031  counter_CAD_2 = counter_CAD_2 + delta_t * 360. * speed / 60.;
4032  } while(counter_CAD_2 <= 180);
4033  complete_post_size--;
4034 
4035  complete_size = size + complete_prev_size + complete_post_size;
4036 
4037  complete_p_cyl = (double*) malloc(complete_size * sizeof(double));
4038  complete_V_cyl = (double*) malloc(complete_size * sizeof(double));
4039  complete_deform = (double*) malloc(complete_size * sizeof(double));
4040  complete_CAD = (double*) malloc(complete_size * sizeof(double));
4041 
4042  CALCULUS_OF_IMP_HP(complete_p_cyl, complete_CAD, p_cyl, V_cyl, complete_V_cyl, complete_deform, &WI_HP, &IMP_HP,
4043  complete_size, complete_prev_size, delta_t, speed, size, IVC, EVO, VTDC,
4044  Cylinder_capacity, PI, Piston_D, S, Crank_L, Connecting_Rod_L, E, Piston_Axis_D, Piston_Crown_H, M_Connecting_Rod,
4045  M_P_R_PA, C_ESteel, C_Mech_Defor, inlet_pres, exhaust_pres);
4046 
4047  CALCULUS_OF_MEAN_VARIABLES(p_cyl, T_cyl, dp_da_cyl, CAD, &pmax, &Tmax, &dp_da_max, &p_exit, &T_exit, size);
4048 
4049  // indicated average power and indicated efficiency
4050  mean_var_exit[0] = IMP_HP * 1e-5;
4051  mean_var_exit[1] = IMP_HP * Cylinder_capacity / (mfuel * 1e-6 * HP);
4052 
4053  // mean variables in.cylinder laws results
4054  mean_var_exit[2] = pmax * 1e-5;
4055  mean_var_exit[3] = Tmax;
4056  mean_var_exit[4] = p_exit * 1e-5;
4057  mean_var_exit[5] = T_exit;
4058  mean_var_exit[6] = dp_da_max * 1e-5;
4059 
4060  // NOx concentration exit
4061  mean_var_exit[7] = mNOx_bowl[size - 3] * 1e6 / (1.587 * (mairIVC + mfuel * 1.e-3));
4062 
4063  mean_var_exit[8] = mO2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4064 
4065  // Species concentration at EVO
4066  (*dataOUT).species_EVO[0] = mN2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4067  (*dataOUT).species_EVO[1] = mO2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4068  (*dataOUT).species_EVO[2] = mCO2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4069  (*dataOUT).species_EVO[3] = mH2O_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4070  (*dataOUT).species_EVO[4] = mNOx_bowl[size - 3] / (mairIVC + mfuel * 1.e-3) / 1.587;
4071  (*dataOUT).species_EVO[5] = mCO_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4072  (*dataOUT).species_EVO[6] = mSOOT_bowl_A[size - 3] / (mairIVC + mfuel * 1.e-3);
4073  (*dataOUT).species_EVO[7] = mHC_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4074  // The mass fractions are adjusted so as to ensure that the summation is 1
4075  Tot = 0.;
4076  for(i = 0; i < 8; i++)
4077  Tot += (*dataOUT).species_EVO[i];
4078  for(i = 0; i < 8; i++)
4079  (*dataOUT).species_EVO[i] /= Tot;
4080 
4081  SOOT_EVO_A = mSOOT_bowl_A[size - 1] * 1.e6 / (mairIVC + mfuel * 1.e-3);
4082 
4083  // Pasamos el valor a FSN
4084  mean_var_exit[10] = YSoot_to_FSN(SOOT_EVO_A);
4085 
4086  SOOT_EVO_C = (evol_mSoot[size - 1] + mSootCil[size - 1]) * 1.e6 / (mairIVC + mfuel * 1.e-3);
4087 
4088  // Pasamos el valor a FSN
4089  mean_var_exit[9] = YSoot_to_FSN(SOOT_EVO_C);
4090  // Anadimos la correccion emperica:
4091  mean_var_exit[9] = 0.0930816 * pow(mean_var_exit[9], 2.68);
4092  // Limitamos el valor a 10
4093  if(mean_var_exit[9] > 10.)
4094  mean_var_exit[9] = 10.;
4095 
4096  heat_transfer[0] = 0;
4097  heat_transfer[1] = 0;
4098  heat_transfer[2] = 0;
4099  heat_transfer[3] = 0;
4100 
4101  for(counter = 0; counter < size - 1; counter++) {
4102  heat_transfer[0] = heat_transfer[0] + min(H_cooler[counter],
4103  H_cooler[counter + 1]) + pow(pow((H_cooler[counter + 1] - H_cooler[counter]) / 2, 2), 0.5);
4104  heat_transfer[1] = heat_transfer[1] + min(Qcylhead[counter],
4105  Qcylhead[counter + 1]) + pow(pow((Qcylhead[counter + 1] - Qcylhead[counter]) / 2, 2), 0.5);
4106  heat_transfer[2] = heat_transfer[2] + min(Qcyl[counter],
4107  Qcyl[counter + 1]) + pow(pow((Qcyl[counter + 1] - Qcyl[counter]) / 2, 2), 0.5);
4108  heat_transfer[3] = heat_transfer[3] + min(Qpis[counter],
4109  Qpis[counter + 1]) + pow(pow((Qpis[counter + 1] - Qpis[counter]) / 2, 2), 0.5);
4110  }
4111 
4112  // instantaneous law cycle control
4113  /*
4114  if((fich=fopen("c:\\ACT_model\\output\\seguimiento_paq_50_3.dat","w"))==NULL){
4115  printf("The element file results could not be opened");
4116  exit(1);
4117  }
4118  strcpy(title,"CAD[deg],time[s],mtotal,mfuel,mfuel_real,mO2,mO2_real,TSD,Tadib");
4119  fprintf(fich, "%s", title);
4120 
4121  for(counter=0;counter<size;counter++){
4122  fprintf(fich, "\n%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf",CAD[counter],time_vector[counter],elementcontrol[7][counter].mtotal,elementcontrol[7][counter].mfuel,elementcontrol[7][counter].mfuel_real,elementcontrol[7][counter].mO2,elementcontrol[7][counter].mO2_real,elementcontrol[7][counter].TSD,elementcontrol[7][counter].Tadib);
4123  }
4124 
4125  fclose(fich);
4126 
4127 
4128  if((fich=fopen("c:\\ACT_model\\output\\Instantaneous_law.dat","w"))==NULL){
4129  printf("The element file results could not be opened");
4130  exit(1);
4131  }
4132  strcpy(title,"CAD[deg],time[s],Pcyl[bar],dpdacyl[bar/deg],HR[J],ROHR[J/s],Tcyl[K],Vcyl[m3],ro_air[kg/m3],uo*[m/s],mtotal_cyl[g],injection[kg/s],accu_mfuel[kg],mO2bowl[g],mO2Vd[g],YO2[-],Tsq[K],Tadib_cyl[K],mNOx[g],XLO[m]");
4133  fprintf(fich, "%s", title);
4134 
4135  for(counter=0;counter<size;counter++){
4136 
4137  fprintf(fich, "\n%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf",CAD[counter],time_vector[counter],(p_cyl[counter]*1e-5),dp_da_cyl[counter],(HRF[counter]*mfuel*42.92),ROHR[counter],\
4138  T_cyl[counter],V_cyl[counter],roair[counter],virtual_velocity[counter],(mtotal_bowl[counter]+mtotal_Vc[counter]),dmf[counter],(acu_dmf[counter]*mfuel*42.92),mO2_bowl[counter],mO2_Vc[counter],((mO2_bowl[counter]+mO2_Vc[counter])/(mtotal_bowl[counter]+mtotal_Vc[counter])),Tsq_cyl[counter],Tadib_cyl[counter],\
4139  (mNOx_bowl[counter]+mNOx_Vc[counter]),XLO[counter]);
4140  }
4141  fclose(fich);
4142 
4143 
4144  if((fich=fopen("c:\\ACT_model\\output\\mtotal.dat","w"))==NULL){
4145  printf("The element file results could not be opened");
4146  exit(1);
4147  }
4148  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4149  fprintf(fich, "%s", title);
4150 
4151  for(counter=0;counter<size;counter++){
4152  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4153  for(i_aux=0;i_aux<15;i_aux++){
4154  fprintf(fich, "%lf,",elementcontrol[i_aux][counter].mtotal);
4155  }
4156  }
4157  fclose(fich);
4158 
4159  if((fich=fopen("c:\\ACT_model\\output\\mO2.dat","w"))==NULL){
4160  printf("The element file results could not be opened");
4161  exit(1);
4162  }
4163  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4164  fprintf(fich, "%s", title);
4165 
4166  for(counter=0;counter<size;counter++){
4167  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4168  for(i_aux=0;i_aux<15;i_aux++){
4169  fprintf(fich, "%lf,",elementcontrol[i_aux][counter].mO2);
4170  }
4171  }
4172  fclose(fich);
4173 
4174  if((fich=fopen("c:\\ACT_model\\output\\mO2_real.dat","w"))==NULL){
4175  printf("The element file results could not be opened");
4176  exit(1);
4177  }
4178  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4179  fprintf(fich, "%s", title);
4180 
4181  for(counter=0;counter<size;counter++){
4182  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4183  for(i_aux=0;i_aux<15;i_aux++){
4184  fprintf(fich, "%lf,",elementcontrol[i_aux][counter].mO2_real);
4185  }
4186  }
4187  fclose(fich);
4188 
4189 
4190 
4191  if((fich=fopen("c:\\ACT_model\\output\\mfuel.dat","w"))==NULL){
4192  printf("The element file results could not be opened");
4193  exit(1);
4194  }
4195  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4196  fprintf(fich, "%s", title);
4197 
4198  for(counter=0;counter<size;counter++){
4199  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4200  for(i_aux=0;i_aux<15;i_aux++){
4201  fprintf(fich, "%lf,",elementcontrol[i_aux][counter].mfuel);
4202  }
4203  }
4204  fclose(fich);
4205 
4206  if((fich=fopen("c:\\ACT_model\\output\\mfuel_real.dat","w"))==NULL){
4207  printf("The element file results could not be opened");
4208  exit(1);
4209  }
4210  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4211  fprintf(fich, "%s", title);
4212 
4213  for(counter=0;counter<size;counter++){
4214  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4215  for(i_aux=0;i_aux<15;i_aux++){
4216  fprintf(fich, "%lf,",elementcontrol[i_aux][counter].mfuel_real);
4217  }
4218  }
4219  fclose(fich);
4220 
4221  if((fich=fopen("c:\\ACT_model\\output\\FI.dat","w"))==NULL){
4222  printf("The element file results could not be opened");
4223  exit(1);
4224  }
4225  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4226  fprintf(fich, "%s", title);
4227 
4228  for(counter=0;counter<size;counter++){
4229  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4230  for(i_aux=0;i_aux<15;i_aux++){
4231 
4232  fprintf(fich, "%lf,",(elementcontrol[i_aux][counter].FI));
4233 
4234  }
4235  }
4236  fclose(fich);
4237 
4238 
4239  if((fich=fopen("c:\\ACT_model\\output\\temperature_SD.dat","w"))==NULL){
4240  printf("The element file results could not be opened");
4241  exit(1);
4242  }
4243  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4244  fprintf(fich, "%s", title);
4245 
4246  for(counter=0;counter<size;counter++){
4247  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4248  for(i_aux=0;i_aux<15;i_aux++){
4249  fprintf(fich, "%lf,",(elementcontrol[i_aux][counter].TSD));
4250  }
4251  }
4252  fclose(fich);
4253  if((fich=fopen("c:\\ACT_model\\output\\temperature_adib.dat","w"))==NULL){
4254  printf("The element file results could not be opened");
4255  exit(1);
4256  }
4257  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4258  fprintf(fich, "%s", title);
4259 
4260  for(counter=0;counter<size;counter++){
4261  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4262  for(i_aux=0;i_aux<15;i_aux++){
4263  fprintf(fich, "%lf,",(elementcontrol[i_aux][counter].Tadib));
4264  }
4265  }
4266  fclose(fich);
4267 
4268 
4269  if((fich=fopen("c:\\ACT_model\\output\\penetration.dat","w"))==NULL){
4270  printf("The element file results could not be opened");
4271  exit(1);
4272  }
4273  strcpy(title,"CAD[deg],time[s],paq_10_1,paq_10_3,paq_10_5,paq_30_1,paq_30_3,paq_30_5,paq_50_1,paq_50_3,paq_50_5,paq_70_1,paq_70_3,paq_70_5,paq_90_1,paq_90_3,paq_90_5");
4274  fprintf(fich, "%s", title);
4275 
4276  for(counter=0;counter<size;counter++){
4277  fprintf(fich, "\n%lf,%lf,",CAD[counter],time_vector[counter]);
4278  for(i_aux=0;i_aux<15;i_aux++){
4279  fprintf(fich, "%lf,",(elementcontrol[i_aux][counter].X));
4280  }
4281  }
4282  fclose(fich);
4283  */
4284 
4285  // exit vectors construction with constant crank angle degree increment
4286  time_vector_exit = (double*) malloc(CAI * sizeof(double));
4287  time_vector_exit[0] = -180 * 60 / (360. * speed);
4288  CAD_exit[0] = -180;
4289 
4290  for(counter = 1; counter < CAI; counter++) {
4291  time_vector_exit[counter] = time_vector_exit[counter - 1] + 60 / ((CAI - 1) * speed);
4292  CAD_exit[counter] = CAD_exit[counter - 1] + 360. / (CAI - 1);
4293  }
4294  // pcyl
4295  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4296  for(counter = 0; counter < size; counter++) {
4297  vector_to_interpolate[counter] = p_cyl[counter];
4298  }
4299  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4300  free(vector_to_interpolate);
4301  for(counter = 0; counter < CAI; counter++) {
4302  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4303  aux_vector[counter] = 0;
4304  }
4305  p_cyl_exit[counter] = aux_vector[counter];
4306  }
4307 
4308  // dp_da_exit
4309  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4310  for(counter = 0; counter < size; counter++) {
4311  vector_to_interpolate[counter] = dp_da_cyl[counter];
4312  // vector_to_interpolate[counter]=PEN_min[counter]*1e5 = 0.; // El factor es para evitar cambio de unidad. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4313  }
4314  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4315  free(vector_to_interpolate);
4316  for(counter = 0; counter < CAI; counter++) {
4317  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4318  aux_vector[counter] = 0;
4319  }
4320  dp_da_cyl_exit[counter] = aux_vector[counter];
4321  }
4322 
4323  // T_cyl
4324  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4325  for(counter = 0; counter < size; counter++) {
4326  vector_to_interpolate[counter] = T_cyl[counter];
4327  }
4328  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4329  free(vector_to_interpolate);
4330  for(counter = 0; counter < CAI; counter++) {
4331  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4332  aux_vector[counter] = 0;
4333  }
4334  T_cyl_exit[counter] = aux_vector[counter];
4335  }
4336 
4337  // HRF
4338  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4339  for(counter = 0; counter < size; counter++) {
4340  vector_to_interpolate[counter] = HRF[counter];
4341  }
4342  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4343  free(vector_to_interpolate);
4344  for(counter = 0; counter < CAI; counter++) {
4345  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4346  aux_vector[counter] = 0;
4347  }
4348  HRF_exit[counter] = aux_vector[counter];
4349  }
4350 
4352  if(aux_mfuel == 0.) {
4353  for(counter = 0; counter < CAI; counter++) {
4354  HRF_exit[counter] = 0.;
4355  }
4356  }
4357 
4358  // HRF_PMX
4359  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4360  for(counter = 0; counter < size; counter++) {
4361  vector_to_interpolate[counter] = HRF_PMX[counter];
4362  }
4363  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4364  free(vector_to_interpolate);
4365  for(counter = 0; counter < CAI; counter++) {
4366  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4367  aux_vector[counter] = 0;
4368  }
4369  (*dataOUT).HRF_PMX[counter] = aux_vector[counter];
4370  }
4371 
4372  // ROHR
4373  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4374  for(counter = 0; counter < size; counter++) {
4375  vector_to_interpolate[counter] = ROHR[counter];
4376  }
4377  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4378  free(vector_to_interpolate);
4379  for(counter = 0; counter < CAI; counter++) {
4380  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4381  aux_vector[counter] = 0;
4382  }
4383  ROHR_exit[counter] = aux_vector[counter];
4384  }
4385 
4386  // H_Cooler
4387  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4388  for(counter = 0; counter < size; counter++) {
4389  vector_to_interpolate[counter] = H_cooler[counter];
4390  }
4391  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4392  free(vector_to_interpolate);
4393  for(counter = 0; counter < CAI; counter++) {
4394  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4395  aux_vector[counter] = 0;
4396  }
4397  H_cooler_exit[counter] = aux_vector[counter] / (delta_t * speed * 360. / 60);
4398  }
4399 
4400  // injection_rate
4401  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4402  for(counter = 0; counter < size; counter++) {
4403  vector_to_interpolate[counter] = dmf[counter];
4404  }
4405  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4406  free(vector_to_interpolate);
4407  for(counter = 0; counter < CAI; counter++) {
4408  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4409  aux_vector[counter] = 0;
4410  }
4411  injection_rate_exit[counter] = aux_vector[counter];
4412  }
4413 
4414  // Accum_injection_rate
4415  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4416  for(counter = 0; counter < size; counter++) {
4417  vector_to_interpolate[counter] = acu_dmf[counter];
4418  // vector_to_interpolate[counter]=PEN_max[counter]/(test_variables[4]*42.92) = 0.; // El factor es para evitar cambio de unidad. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4419  }
4420  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4421  free(vector_to_interpolate);
4422  for(counter = 0; counter < CAI; counter++) {
4423  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4424  aux_vector[counter] = 0;
4425  }
4426  accum_injection_rate_exit[counter] = aux_vector[counter];
4427  }
4428 
4429  // evol_Soot
4430  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4431  for(counter = 0; counter < size; counter++) {
4432  vector_to_interpolate[counter] = evol_mSoot[counter];
4433  }
4434  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4435  free(vector_to_interpolate);
4436  for(counter = 0; counter < CAI; counter++) {
4437  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4438  aux_vector[counter] = 0;
4439  }
4440  (*dataOUT).evol_Soot[counter] = aux_vector[counter];
4441  }
4442 
4443  // evol_Soot_CIL
4444  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4445  for(counter = 0; counter < size; counter++) {
4446  vector_to_interpolate[counter] = mSootCil[counter];
4447  }
4448  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4449  free(vector_to_interpolate);
4450  for(counter = 0; counter < CAI; counter++) {
4451  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4452  aux_vector[counter] = 0;
4453  }
4454  (*dataOUT).evol_Soot_CIL[counter] = aux_vector[counter];
4455  }
4456 
4457  // evol_Radiacion
4458  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4459  for(counter = 0; counter < size; counter++) {
4460  vector_to_interpolate[counter] = evol_Radiacion[counter];
4461  }
4462  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4463  free(vector_to_interpolate);
4464  for(counter = 0; counter < CAI; counter++) {
4465  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4466  aux_vector[counter] = 0;
4467  }
4468  (*dataOUT).evol_Radiacion[counter] = aux_vector[counter];
4469  }
4470 
4471  // evol_LOL
4472  vector_to_interpolate = (double*) malloc(size * sizeof(double));
4473  for(counter = 0; counter < size; counter++) {
4474  vector_to_interpolate[counter] = XLO[counter];
4475  }
4476  FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4477  free(vector_to_interpolate);
4478  for(counter = 0; counter < CAI; counter++) {
4479  if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4480  aux_vector[counter] = 0;
4481  }
4482  (*dataOUT).evol_LOL[counter] = aux_vector[counter];
4483  }
4484 
4487  free(CAD);
4488  free(time_vector);
4489  free(EOI_IM);
4490  free(SOI_IM);
4491  free(SOC_IM);
4492 
4493  free(mtotal_bowl);
4494  free(mO2_bowl);
4495  free(mN2_bowl);
4496  free(mCO2_bowl);
4497  free(mH2O_bowl);
4498  free(mHC_bowl);
4499  free(Tsq_cyl);
4500  free(Tadib_cyl);
4501  free(mNOx_bowl);
4502  free(mCO_bowl);
4503  free(HC_bowl);
4504  free(mSOOT_bowl_A);
4505  free(mSOOT_bowl_B);
4506  free(mSOOT_bowl_C);
4507 
4508  free(mtotal_Vc);
4509  free(mO2_Vc);
4510  free(mN2_Vc);
4511  free(mCO2_Vc);
4512  free(mH2O_Vc);
4513  free(mNOx_Vc);
4514  free(mSOOT_Vc_A);
4515  free(mSOOT_Vc_B);
4516  free(mSOOT_Vc_C);
4517  free(mHC_Vc);
4518  free(dmtotal_Gfactor);
4519 
4520  free(evol_mSoot);
4521  free(mSootCil);
4522  free(evol_Radiacion);
4523 
4524  free(dmf);
4525  free(acu_dmf);
4526 
4527  free(virtual_velocity);
4528  free(inj_velocity);
4529 
4530  free(HRF);
4531  free(HRF_PMX);
4532  free(ROHR);
4533 
4534  free(T_cyl);
4535  free(acu_Mbb);
4536  free(V_cyl);
4537  free(Mbb);
4538  free(Yair);
4539  free(Yburned);
4540  free(Yfuel);
4541  free(U);
4542  free(CV);
4543  free(H_cooler);
4544  free(Qcylhead);
4545  free(Qcyl);
4546  free(Qpis);
4547  free(H);
4548  free(defor);
4549  free(Rmixture);
4550  free(roair);
4551  free(p_cyl);
4552  free(dp_da_cyl);
4553  free(aux_vector);
4554  free(XLO);
4555  free(PEN_min);
4556  free(PEN_max);
4557 
4558  for(counter = 0; counter < 52; counter++)
4559  free(YNOeq[counter]);
4560  free(YNOeq);
4561  for(counter = 0; counter < 52; counter++)
4562  free(KdYNO[counter]);
4563  free(KdYNO);
4564 
4565  free(num_i_IM);
4566 
4567  for(counter = 0; counter < inj_num; counter++)
4568  free(POI_IM[counter]);
4569  free(POI_IM);
4570  for(counter = 0; counter < inj_num; counter++)
4571  free(POC_IM[counter]);
4572  free(POC_IM);
4573  for(counter = 0; counter < inj_num; counter++)
4574  free(mfuel_ij_IM[counter]);
4575  free(mfuel_ij_IM);
4576  for(counter = 0; counter < inj_num; counter++)
4577  free(mfuel_i_IM[counter]);
4578  free(mfuel_i_IM);
4579 
4580  for(counter = 0; counter < 15; counter++)
4581  free(elementcontrol[counter]);
4582  free(elementcontrol);
4583 
4584  free(mixture_correction);
4585 
4586  for(counter = 0; counter < inj_num; counter++)
4587  free(element[counter]);
4588  free(element);
4589  for(counter = 0; counter < inj_num; counter++)
4590  free(realelement[counter]);
4591  free(realelement);
4592 
4593  free(complete_p_cyl);
4594  free(complete_V_cyl);
4595  free(complete_deform);
4596  free(complete_CAD);
4597 
4598  free(time_vector_exit);
4599 
4600 }
4601 // ------------------------------------------------------
4602 
4603 double YSoot_to_FSN(double YSoot) {
4604  double x1 = 0., x2 = 0., x = 0., y = 0.;
4605  // Busca el FSN implentando un metodo de la biseccion.
4606  if(YSoot < 0.)
4607  x = 0.;
4608  else {
4609  if(YSoot > 4550.)
4610  x = 10.;
4611  else {
4612  x1 = 0;
4613  x2 = 10;
4614  while((x2 - x1) > 0.01) {
4615  x = (x1 + x2) / 2.;
4616  y = 1 / 1.2 / 0.405 * 4.95 * x * exp(0.38 * x);
4617  if(y > YSoot) {
4618  x2 = x;
4619  } else {
4620  x1 = x;
4621  }
4622  }
4623  }
4624  }
4625  return (x);
4626 }
stControlElementComposition
Definition: ACT_Sub_DLL.h:49
Constantes.h
stRadArray
Definition: ACT_Sub_DLL.h:35
sOUTtype
Definition: ACT_Sub_DLL.h:26
sINtype
Definition: ACT_Sub_DLL.h:22