cloudy
trunk
Main Page
Related Pages
Data Structures
Files
File List
Globals
All
Data Structures
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
source
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"
,
307
radius
.
depth_mid_zone
,
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
1.8.1.2