cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
heat_punch.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*HeatPunch punch contributors to local heating, with punch heat command, called by punch_do */
4 #include "cddefines.h"
5 #include "thermal.h"
6 #include "radius.h"
7 #include "conv.h"
8 #include "lines_service.h"
9 #include "dense.h"
10 #include "taulines.h"
11 #include "phycon.h"
12 #include "elementnames.h"
13 #include "dynamics.h"
14 #include "punch.h"
15 
16 /* limit for number of heat agents that are saved */
17 /* limit to number to print */
18 static const int IPRINT= 9;
19 
20 /*HeatPunch punch contributors to local heating, with punch heat command,
21  * called by punch_do */
22 void HeatPunch(FILE* io)
23 {
24  char **chLabel,
25  chLbl[11];
26  bool lgHeatLine;
27  int nFail;
28  long int i,
29  ipStrong,
30  ipnt,
31  *ipOrdered,
32  *ipsave,
33  j,
34  *jpsave,
35  k,
36  level;
37  double CS,
38  ColHeat,
39  EscP,
40  Pump,
41  Strong,
42  TauIn,
43  cool_total,
44  heat_total;
45  realnum *save;
46 
47  DEBUG_ENTRY( "HeatPunch()" );
48 
49  save = (realnum *) CALLOC(LIMELM*LIMELM, sizeof(realnum));
50  ipsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
51  jpsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
52  ipOrdered = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
53  chLabel = (char **) CALLOC(LIMELM*LIMELM, sizeof(char *));
54 
55  for( i=0; i<LIMELM*LIMELM; ++i )
56  {
57  ipsave[i] = INT_MIN;
58  jpsave[i] = INT_MIN;
59  save[i] = -FLT_MAX;
60  chLabel[i] = (char *) CALLOC(10, sizeof(char));
61  }
62 
63  cool_total = thermal.ctot;
64  heat_total = thermal.htot;
65 
66  /* >>chng 06 mar 17, comment out following block and replace with this
67  * removing dynamics heating & cooling and report only physical
68  * heating and cooling
69  * NB the heating and cooling as punched no longer need be
70  * equal for a converged model */
71  cool_total -= dynamics.Cool;
72  heat_total -= dynamics.Heat;
73 # if 0
74  if(dynamics.Cool > dynamics.Heat)
75  {
76  cool_total -= dynamics.Heat;
77  heat_total -= dynamics.Heat;
78  }
79  else
80  {
81  cool_total -= dynamics.Cool;
82  heat_total -= dynamics.Cool;
83  }
84 # endif
85 
86  ipnt = 0;
87 
88  /* heat sources are saved in a 2d square array */
89  /* WeakHeatCool set with 'set weakheatcool' command
90  * default is 0.05 */
91  for( i=0; i < LIMELM; i++ )
92  {
93  for( j=0; j < LIMELM; j++ )
94  {
95  if( thermal.heating[i][j]/heat_total > SMALLFLOAT )
96  {
97  ipsave[ipnt] = i;
98  jpsave[ipnt] = j;
99  save[ipnt] = (realnum)(thermal.heating[i][j]/heat_total);
100  ipnt++;
101  }
102  }
103  }
104 
105  /* now check for possible line heating not counted in 1,23
106  * there should not be any significant heating source here
107  * since they would not be counted in derivative correctly */
108  for( i=0; i < thermal.ncltot; i++ )
109  {
110  if( thermal.heatnt[i]/heat_total > punch.WeakHeatCool )
111  {
112  realnum awl;
113  awl = thermal.collam[i];
114  /* force to save wavelength convention as printout */
115  if( awl > 100000 )
116  awl /= 10000;
117  fprintf( io, " Negative coolant was %s %.2f %.2e\n",
118  thermal.chClntLab[i], awl, thermal.heatnt[i]/heat_total );
119  }
120  }
121 
122  if( !conv.lgConvTemp )
123  {
124  fprintf( io, "#>>>> Temperature not converged.\n" );
125  }
126  else if( !conv.lgConvEden )
127  {
128  fprintf( io, "#>>>> Electron density not converged.\n" );
129  }
130  else if( !conv.lgConvIoniz )
131  {
132  fprintf( io, "#>>>> Ionization not converged.\n" );
133  }
134  else if( !conv.lgConvPres )
135  {
136  fprintf( io, "#>>>> Pressure not converged.\n" );
137  }
138 
139  /* this is mainly to keep the compiler from flagging possible paths
140  * with j not being set */
141  i = INT_MIN;
142  j = INT_MIN;
143  /* following loop tries to identify strongest agents and turn to labels */
144  for( k=0; k < ipnt; k++ )
145  {
146  /* generate labels that make sense in printout
147  * if not identified with a specific agent, will print indices as [i][j],
148  * heating is thermal.heating[i][j] */
149  i = ipsave[k];
150  j = jpsave[k];
151  /* i >= j indicates agent is one of the first LIMELM elements */
152  if( i >= j )
153  {
154  if( dense.xIonDense[i][j] == 0. && thermal.heating[i][j]>SMALLFLOAT )
155  fprintf(ioQQQ,"DISASTER assert about to be thrown - search for hit it\n");
156  /*fprintf(ioQQQ,"DEBUG %li %li %.2e %.2e\n", i , j ,
157  dense.xIonDense[i][j],
158  thermal.heating[i][j]);*/
159  ASSERT( dense.xIonDense[i][j] > 0. || thermal.heating[i][j]<SMALLFLOAT );
160  /* this is case of photoionization of atom or ion */
161  strcpy( chLabel[k], elementnames.chElementSym[i] );
162  strcat( chLabel[k], elementnames.chIonStage[j] );
163  }
164  /* notice that in test i and j are swapped from order in heating array */
165  else if( i == 0 && j == 1 )
166  {
167  /* photoionization from all excited states of Hydrogenic species,
168  * heating[0][1] */
169  strcpy( chLabel[k], "Hn=2" );
170  }
171  else if( i == 0 && j == 3 )
172  {
173  /* collisional ionization of all iso-seq from all levels,
174  * heating[0][3] */
175  strcpy( chLabel[k], "Hion" );
176  }
177  else if( i == 0 && j == 7 )
178  {
179  /* UTA ionization heating[0][7] */
180  strcpy( chLabel[k], " UTA" );
181  }
182  else if( i == 0 && j == 8 )
183  {
184  /* thermal.heating[0][8] is heating due to collisions within
185  * X of H2, code var hmi.HeatH2Dexc_used, hmi.HeatH2Dexc_BigH2,
186  * hmi.HeatH2Dexc_TH85 */
187  strcpy( chLabel[k], "H2vH" );
188  }
189  else if( i == 0 && j == 17 )
190  {
191  /* thermal.heating[0][17] is heating due to photodissociation
192  * heating of X within H2,
193  * code var hmi.HeatH2Dish_used */
194  strcpy( chLabel[k], "H2dH" );
195  }
196  else if( i == 0 && j == 9 )
197  {
198  /* CO dissociation, co.CODissHeat, heating[0][9] */
199  strcpy( chLabel[k], "COds" );
200  }
201  else if( i == 0 && j == 20 )
202  {
203  /* extra heat thermal.heating[0][20]*/
204  strcpy( chLabel[k], "extH" );
205  }
206  else if( i == 0 && j == 11 )
207  {
208  /* free free heating, heating[0][11] */
209  strcpy( chLabel[k], "H FF" );
210  }
211  else if( i == 0 && j == 12 )
212  {
213  /* heating line, that was supposed to cool, heating[0][12] */
214  strcpy( chLabel[k], "Clin" );
215  }
216  else if( i == 0 && j == 13 )
217  {
218  /* grain photoionization, heating[0][13] */
219  strcpy( chLabel[k], "GrnP" );
220  }
221  else if( i == 0 && j == 14 )
222  {
223  /* grain collisions, heating[0][14] */
224  strcpy( chLabel[k], "GrnC" );
225  }
226  else if( i == 0 && j == 15 )
227  {
228  /* H- heating, heating[0][15] */
229  strcpy( chLabel[k], "H- " );
230  }
231  else if( i == 0 && j == 16 )
232  {
233  /* H2+ heating, heating[0][16] */
234  strcpy( chLabel[k], "H2+ " );
235  }
236  else if( i == 0 && j == 18 )
237  {
238  /* H2 photoionization heating, heating[0][18] */
239  strcpy( chLabel[k], "H2ph" );
240  }
241  else if( i == 0 && j == 19 )
242  {
243  /* Compton heating, heating[0][19] */
244  strcpy( chLabel[k], "Comp" );
245  }
246  else if( i == 0 && j == 22 )
247  {
248  /* line heating, heating[0][22] */
249  strcpy( chLabel[k], "line" );
250  }
251  else if( i == 0 && j == 23 )
252  {
253  /* iso-sequence line heating - all elements together,
254  * heating[0][23] */
255  strcpy( chLabel[k], "Hlin" );
256  }
257  else if( i == 0 && j == 24 )
258  {
259  /* charge transfer heating, heating[0][24] */
260  strcpy( chLabel[k], "ChaT" );
261  }
262  else if( i == 1 && j == 3 )
263  {
264  /* helium triplet line heating, heating[1][3] */
265  strcpy( chLabel[k], "He3l" );
266  }
267  else if( i == 1 && j == 5 )
268  {
269  /* advective heating, heating[1][5] */
270  strcpy( chLabel[k], "adve" );
271  }
272  else if( i == 1 && j == 6 )
273  {
274  /* cosmic ray heating thermal.heating[1][6]*/
275  strcpy( chLabel[k], "CR H" );
276  }
277  else if( i == 25 && j == 27 )
278  {
279  /* Fe 2 line heating, heating[25][27] */
280  strcpy( chLabel[k], "Fe 2" );
281  }
282  else
283  {
284  sprintf( chLabel[k], "[%ld][%ld]" , i , j );
285  }
286  }
287 
288  /* now sort by decreasing importance */
289  /*spsort netlib routine to sort array returning sorted indices */
290  spsort(
291  /* input array to be sorted */
292  save,
293  /* number of values in x */
294  ipnt,
295  /* permutation output array */
296  ipOrdered,
297  /* flag saying what to do - 1 sorts into increasing order, not changing
298  * the original routine */
299  -1,
300  /* error condition, should be 0 */
301  &nFail);
302 
303  /*>>chng 06 jun 06, change start of punch to give same info as cooling
304  * as per comment by Yumihiko Tsuzuki */
305  /* begin the print out with zone number, total heating and cooling */
306  fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e",
308  phycon.te,
309  heat_total,
310  cool_total );
311 
312  /* we only want to print the IPRINT strongest of the group */
313  ipnt = MIN2( ipnt , IPRINT );
314 
315  for( k=0; k < ipnt; k++ )
316  {
317  int ip = ipOrdered[k];
318  i = ipsave[ip];
319  j = jpsave[ip];
320  ASSERT( i<LIMELM && j<LIMELM );
321  if(k > 4 && thermal.heating[i][j]/heat_total < punch.WeakHeatCool)
322  break;
323  fprintf( io, "\t%s\t%.7f ",
324  chLabel[ip], save[ip] );
325  }
326  fprintf( io, " \n" );
327 
328  /* a negative pointer in the heating array is probably a problem,
329  * indicating that some line has become a heat source */
330  lgHeatLine = false;
331 
332  /* check if any lines were major heat sources */
333  for( i=0; i < ipnt; i++ )
334  {
335  /* heating[22][0] is line heating - identify line if important */
336  if( ipsave[i] == 0 && jpsave[i] == 22 )
337  lgHeatLine = true;
338  }
339 
340  if( lgHeatLine )
341  {
342  /* a line was a major heat source - identify it */
343  FndLineHt(&level,&ipStrong,&Strong);
344  if( Strong/heat_total > 0.005 )
345  {
346  if( level == 1 )
347  {
348  strcpy( chLbl, chLineLbl(&TauLines[ipStrong] ) );
349  TauIn = TauLines[ipStrong].Emis->TauIn;
350  Pump = TauLines[ipStrong].Emis->pump;
351  EscP = TauLines[ipStrong].Emis->Pesc;
352  CS = TauLines[ipStrong].Coll.cs;
353  /* ratio of line to total heating */
354  ColHeat = TauLines[ipStrong].Coll.heat/
355  heat_total;
356  }
357  else if( level == 2 )
358  {
359  strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]) );
360  TauIn = TauLine2[ipStrong].Emis->TauIn;
361  Pump = TauLine2[ipStrong].Emis->pump;
362  EscP = TauLine2[ipStrong].Emis->Pesc;
363  CS = TauLine2[ipStrong].Coll.cs;
364  ColHeat = TauLine2[ipStrong].Coll.heat/
365  heat_total;
366  }
367  else if( level == 3 )
368  /* C12O16 */
369  {
370  strcpy( chLbl, chLineLbl(&C12O16Rotate[ipStrong]) );
371  TauIn = C12O16Rotate[ipStrong].Emis->TauIn;
372  Pump = C12O16Rotate[ipStrong].Emis->pump;
373  EscP = C12O16Rotate[ipStrong].Emis->Pesc;
374  CS = C12O16Rotate[ipStrong].Coll.cs;
375  ColHeat = C12O16Rotate[ipStrong].Coll.heat/
376  heat_total;
377  }
378  else if( level == 4 )
379  /* C13O16 */
380  {
381  strcpy( chLbl, chLineLbl(&C13O16Rotate[ipStrong]) );
382  TauIn = C13O16Rotate[ipStrong].Emis->TauIn;
383  Pump = C13O16Rotate[ipStrong].Emis->pump;
384  EscP = C13O16Rotate[ipStrong].Emis->Pesc;
385  CS = C13O16Rotate[ipStrong].Coll.cs;
386  ColHeat = C13O16Rotate[ipStrong].Coll.heat/
387  heat_total;
388  }
389  else if( level == 5 )
390  /* hyperfine levels */
391  {
392  strcpy( chLbl, chLineLbl(&HFLines[ipStrong]) );
393  TauIn = HFLines[ipStrong].Emis->TauIn;
394  Pump = HFLines[ipStrong].Emis->pump;
395  EscP = HFLines[ipStrong].Emis->Pesc;
396  CS = HFLines[ipStrong].Coll.cs;
397  ColHeat = HFLines[ipStrong].Coll.heat/
398  heat_total;
399  }
400  else
401  {
402  fprintf( ioQQQ, " HeatPunch insane level%4ld\n",
403  level );
404  cdEXIT(EXIT_FAILURE);
405  }
406  fprintf( io, " LHeat lv%2ld %10.10s TIn%10.2e Pmp%9.1e EscP%9.1e CS%9.1e Hlin/tot%10.2e\n",
407  level, chLbl, TauIn, Pump, EscP, CS, ColHeat );
408  }
409  }
410  for( i=0; i<LIMELM*LIMELM; ++i )
411  {
412  free(chLabel[i]);
413  }
414 
415  free(chLabel);
416  free(ipOrdered);
417  free(jpsave);
418  free(ipsave);
419  free(save);
420  return;
421 }

Generated for cloudy by doxygen 1.8.1.1