root/branches/newmole/source/atom_seq_boron.cpp

Revision 2542, 15.9 kB (checked in by rjrw, 3 weeks ago)

Merged from trunk r2448:2541

  • Property svn:eol-style set to native
Line 
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/*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
4#include "cddefines.h"
5#include "cooling.h"
6#include "thermal.h"
7#include "dense.h"
8#include "atoms.h"
9
10void AtomSeqBoron(
11        /* indices for all lines are on the C scale since they will be stuffed into
12         * C arrays.  so, t10 refers to the 2-1 transition */
13        transition * t10,
14        transition * t20,
15        transition * t30,
16        transition * t21,
17        transition * t31,
18        transition * t41,
19        double cs40,
20        double cs32,
21        double cs42,
22        double cs43,
23        /* pump rate s-1 due to UV permitted lines */
24        double pump_rate ,
25        /* string used to identify calling program in case of error */
26        const char *chLabel
27        )
28{
29
30        /* this routine has three possible returns:
31         * abundance is zero
32         * too cool for full 5-level atom, but still do ground term
33         * full solution */
34
35        /* boron sequence is now a five level atom */
36#       define  N_SEQ_BORON     5
37        static double 
38                **AulEscp ,
39                **col_str ,
40                **AulDest,
41                /* AulPump[low][high] is rate (s^-1) from lower to upper level */
42                **AulPump,
43                **CollRate,
44                *pops,
45                *create,
46                *destroy,
47                *depart,
48                /* statistical weight */
49                *stat ,
50                /* excitation energies in kelvin */
51                *excit;
52
53        double b_cooling,
54          dCoolDT;
55        double EnrLU, EnrUL;
56        realnum abundan;
57
58        static bool lgFirst=true,
59                lgZeroPop;
60        int i , j;
61        int 
62                /* flag to signal negative level populations */
63                lgNegPop;
64                /* flag to turn on debug print in atom_levelN */
65        bool lgDeBug;
66
67        DEBUG_ENTRY( "AtomSeqBoron()" );
68
69        if( lgFirst )
70        {
71                /* will never do this again */
72                lgFirst = false;
73                /* allocate the 1D arrays*/
74                excit = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) );
75                stat = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) );
76                pops = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) );
77                create = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) );
78                destroy = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) );
79                depart = (double *)MALLOC( sizeof(double)*(N_SEQ_BORON) );
80                /* create space for the 2D arrays */
81                AulPump = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)));
82                CollRate = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)));
83                AulDest = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)));
84                AulEscp = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)));
85                col_str = ((double **)MALLOC((N_SEQ_BORON)*sizeof(double *)));
86                for( i=0; i<(N_SEQ_BORON); ++i )
87                {
88                        AulPump[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )));
89                        CollRate[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )));
90                        AulDest[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )));
91                        AulEscp[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )));
92                        col_str[i] = ((double *)MALLOC((N_SEQ_BORON)*sizeof(double )));
93                }
94        }
95
96        /* total abundance of this species */
97        abundan = dense.xIonDense[ t10->Hi->nelem -1][t10->Hi->IonStg-1];
98
99        /** \todo       2       use TransitionZero here */
100        if( abundan <= 0. )
101        {
102                /* this branch, no abundance of ion */
103                t10->Lo->Pop = 0.;
104                t20->Lo->Pop = 0.;
105                t30->Lo->Pop = 0.;
106                t21->Lo->Pop = 0.;
107                t31->Lo->Pop = 0.;
108                t41->Lo->Pop = 0.;
109
110                t10->Emis->PopOpc = 0.;
111                t20->Emis->PopOpc = 0.;
112                t30->Emis->PopOpc = 0.;
113                t21->Emis->PopOpc = 0.;
114                t31->Emis->PopOpc = 0.;
115                t41->Emis->PopOpc = 0.;
116
117                t10->Hi->Pop = 0.;
118                t20->Hi->Pop = 0.;
119                t30->Hi->Pop = 0.;
120                t21->Hi->Pop = 0.;
121                t31->Hi->Pop = 0.;
122                t41->Hi->Pop = 0.;
123
124                t10->Emis->xIntensity = 0.;
125                t20->Emis->xIntensity = 0.;
126                t30->Emis->xIntensity = 0.;
127                t21->Emis->xIntensity = 0.;
128                t31->Emis->xIntensity = 0.;
129                t41->Emis->xIntensity = 0.;
130
131                t10->Coll.cool = 0.;
132                t20->Coll.cool = 0.;
133                t30->Coll.cool = 0.;
134                t21->Coll.cool = 0.;
135                t31->Coll.cool = 0.;
136                t41->Coll.cool = 0.;
137
138                t10->Emis->phots = 0.;
139                t20->Emis->phots = 0.;
140                t30->Emis->phots = 0.;
141                t21->Emis->phots = 0.;
142                t31->Emis->phots = 0.;
143                t41->Emis->phots = 0.;
144
145                t10->Emis->ColOvTot = 0.;
146                t20->Emis->ColOvTot = 0.;
147                t30->Emis->ColOvTot = 0.;
148                t21->Emis->ColOvTot = 0.;
149                t31->Emis->ColOvTot = 0.;
150                t41->Emis->ColOvTot = 0.;
151
152                t10->Coll.heat = 0.;
153                t20->Coll.heat = 0.;
154                t30->Coll.heat = 0.;
155                t21->Coll.heat = 0.;
156                t31->Coll.heat = 0.;
157                t41->Coll.heat = 0.;
158
159                CoolAdd( chLabel, t10->WLAng , 0.);
160                CoolAdd( chLabel, t20->WLAng , 0.);
161                CoolAdd( chLabel, t30->WLAng , 0.);
162                CoolAdd( chLabel, t21->WLAng , 0.);
163                CoolAdd( chLabel, t31->WLAng , 0.);
164                CoolAdd( chLabel, t41->WLAng , 0.);
165
166                /* level populations */
167                /* LIMLEVELN is the dimension of the atoms vectors */
168                ASSERT( N_SEQ_BORON <= LIMLEVELN);
169                for( i=0; i < N_SEQ_BORON; i++ )
170                {
171                        atoms.PopLevels[i] = 0.;
172                        atoms.DepLTELevels[i] = 1.;
173                }
174                return;
175        }
176
177        ASSERT( t10->Coll.col_str > 0.);
178        ASSERT( t20->Coll.col_str > 0.);
179        ASSERT( t30->Coll.col_str > 0.);
180        ASSERT( t21->Coll.col_str > 0.);
181        ASSERT( t31->Coll.col_str > 0.);
182        ASSERT( t41->Coll.col_str > 0.);
183        ASSERT( cs40>0.);
184        ASSERT( cs32>0.);
185        ASSERT( cs42>0.);
186        ASSERT( cs43>0.);
187
188        /* all elements are used, and must be set to zero if zero */
189        for( i=0; i < N_SEQ_BORON; i++ )
190        {
191                create[i] = 0.;
192                destroy[i] = 0.;
193                for( j=0; j < N_SEQ_BORON; j++ )
194                {
195                        /*data[j][i] = -1e33;*/
196                        AulEscp[j][i] = 0.;
197                        AulDest[j][i] = 0.;
198                        AulPump[j][i] = 0.;
199                        col_str[j][i] = 0.;
200                }
201        }
202
203        /* statistical weights */
204        stat[0] = t10->Lo->g;
205        stat[1] = t10->Hi->g;
206        stat[2] = t20->Hi->g;
207        stat[3] = t30->Hi->g;
208        stat[4] = t41->Hi->g;
209        ASSERT( stat[0]>0. && stat[1]>0. &&stat[2]>0. &&stat[3]>0. &&stat[4]>0.);
210        ASSERT( fabs(t10->Lo->g/2.-1.) < FLT_EPSILON);
211        ASSERT( fabs(t10->Hi->g/4.-1.) < FLT_EPSILON);
212        ASSERT( fabs(t20->Lo->g/2.-1.) < FLT_EPSILON);
213        ASSERT( fabs(t20->Hi->g/2.-1.) < FLT_EPSILON);
214        ASSERT( fabs(t30->Lo->g/2.-1.) < FLT_EPSILON);
215        ASSERT( fabs(t30->Hi->g/4.-1.) < FLT_EPSILON);
216        ASSERT( fabs(t21->Lo->g/4.-1.) < FLT_EPSILON);
217        ASSERT( fabs(t21->Hi->g/2.-1.) < FLT_EPSILON);
218        ASSERT( fabs(t31->Lo->g/4.-1.) < FLT_EPSILON);
219        ASSERT( fabs(t31->Hi->g/4.-1.) < FLT_EPSILON);
220        ASSERT( fabs(t41->Lo->g/4.-1.) < FLT_EPSILON);
221        ASSERT( fabs(t41->Hi->g/6.-1.) < FLT_EPSILON);
222
223        /* excitation energy of each level relative to ground, in Kelvin */
224        excit[0] = 0.;
225        excit[1] = t10->EnergyK;
226        excit[2] = t20->EnergyK;
227        excit[3] = t30->EnergyK;
228        excit[4] = t41->EnergyK + t10->EnergyK;
229        ASSERT( excit[1]>0. &&excit[2]>0. &&excit[3]>0. &&excit[4]>0.);
230
231        /* fill in Einstein As, collision strengths, pumping rates */
232        AulEscp[1][0] = t10->Emis->Aul*(t10->Emis->Pesc + t10->Emis->Pelec_esc);
233        AulDest[1][0] = t10->Emis->Aul*t10->Emis->Pdest;
234        col_str[1][0] = t10->Coll.col_str;
235        AulPump[0][1] = t10->Emis->pump;
236
237        /* add FUV pump transitions to this pump rate */
238        AulPump[0][1] += pump_rate;
239
240        AulEscp[2][0] = t20->Emis->Aul*(t20->Emis->Pesc + t20->Emis->Pelec_esc);
241        AulDest[2][0] = t20->Emis->Aul*t20->Emis->Pdest;
242        col_str[2][0] = t20->Coll.col_str;
243        AulPump[0][2] = t20->Emis->pump;
244
245        AulEscp[3][0] = t30->Emis->Aul*(t30->Emis->Pesc + t30->Emis->Pelec_esc);
246        AulDest[3][0] = t30->Emis->Aul*t30->Emis->Pdest;
247        col_str[3][0] = t30->Coll.col_str;
248        AulPump[0][3] = t30->Emis->pump;
249
250        AulEscp[4][0] = 1e-8;/* made up trans prob */
251        AulDest[4][0] = 0.;
252        col_str[4][0] = cs40;
253        AulPump[0][4] = 0.;
254
255        AulEscp[2][1] = t21->Emis->Aul*(t21->Emis->Pesc + t21->Emis->Pelec_esc);
256        AulDest[2][1] = t21->Emis->Aul*t21->Emis->Pdest;
257        col_str[2][1] = t21->Coll.col_str;
258        AulPump[1][2] = t21->Emis->pump;
259
260        AulEscp[3][1] = t31->Emis->Aul*(t31->Emis->Pesc + t31->Emis->Pelec_esc);
261        AulDest[3][1] = t31->Emis->Aul*t31->Emis->Pdest;
262        col_str[3][1] = t31->Coll.col_str;
263        AulPump[1][3] = t31->Emis->pump;
264
265        AulEscp[4][1] = t41->Emis->Aul*(t41->Emis->Pesc + t41->Emis->Pelec_esc);
266        AulDest[4][1] = t41->Emis->Aul*t41->Emis->Pdest;
267        col_str[4][1] = t41->Coll.col_str;
268        AulPump[1][4] = t41->Emis->pump;
269
270        AulEscp[3][2] = 1e-8;/* made up trans prob */
271        AulDest[3][2] = 0.;
272        col_str[3][2] = cs32;
273        AulPump[2][3] = 0.;
274
275        AulEscp[4][2] = 1e-8;/* made up trans prob */
276        AulDest[4][2] = 0.;
277        col_str[4][2] = cs42;
278        AulPump[2][4] = 0.;
279
280        AulEscp[4][3] = 1e-8;/* made up trans prob */
281        AulDest[4][3] = 0.;
282        col_str[4][3] = cs43;
283        AulPump[3][4] = 0.;
284
285        lgDeBug = false;
286
287        /* lgNegPop positive if negative pops occurred, negative if too cold */
288        atom_levelN(N_SEQ_BORON,
289                abundan,
290                stat,
291                excit,
292                'K',
293                pops,
294                depart,
295                &AulEscp,
296                &col_str,
297                &AulDest,
298                &AulPump,
299                &CollRate,
300                create,
301                destroy,
302                false,/* say atom_levelN should evaluate coll rates from cs */
303                &b_cooling,
304                &dCoolDT,
305                chLabel,
306                &lgNegPop,
307                &lgZeroPop,
308                lgDeBug );/* option to print stuff - set to true for debug printout */
309
310        /* this returned lgNegPop < 0 in case were gas too cool to excite highest level.
311         * in this case we still want to evaluate the cooling due to the split ground term
312         * in this sequence */
313        if( lgNegPop < 0 )
314        {
315                /* it was too cool for atom_levelN to do whole atom, but let's still do the
316                 * split ground term - this routine adds cooling and its temperature
317                 * derivative to total cooling and dD/dT */
318                atom_level2(t10);
319
320                /* now put in pop for absorption into higher levels, for line optical depths,
321                 * PopLevels was evaluated by atom_level2 */
322                t20->Lo->Pop = atoms.PopLevels[0];
323                t30->Lo->Pop = atoms.PopLevels[0];
324                t21->Lo->Pop = atoms.PopLevels[1];
325                t31->Lo->Pop = atoms.PopLevels[1];
326                t41->Lo->Pop = atoms.PopLevels[1];
327
328                t20->Emis->PopOpc = atoms.PopLevels[0];
329                t30->Emis->PopOpc = atoms.PopLevels[0];
330                t21->Emis->PopOpc = atoms.PopLevels[1];
331                t31->Emis->PopOpc = atoms.PopLevels[1];
332                t41->Emis->PopOpc = atoms.PopLevels[1];
333
334                t20->Hi->Pop = 0.;
335                t30->Hi->Pop = 0.;
336                t21->Hi->Pop = 0.;
337                t31->Hi->Pop = 0.;
338                t41->Hi->Pop = 0.;
339
340                t20->Emis->xIntensity = 0.;
341                t30->Emis->xIntensity = 0.;
342                t21->Emis->xIntensity = 0.;
343                t31->Emis->xIntensity = 0.;
344                t41->Emis->xIntensity = 0.;
345
346                t20->Coll.cool = 0.;
347                t30->Coll.cool = 0.;
348                t21->Coll.cool = 0.;
349                t31->Coll.cool = 0.;
350                t41->Coll.cool = 0.;
351
352                t20->Emis->phots = 0.;
353                t30->Emis->phots = 0.;
354                t21->Emis->phots = 0.;
355                t31->Emis->phots = 0.;
356                t41->Emis->phots = 0.;
357
358                t20->Emis->ColOvTot = 0.;
359                t30->Emis->ColOvTot = 0.;
360                t21->Emis->ColOvTot = 0.;
361                t31->Emis->ColOvTot = 0.;
362                t41->Emis->ColOvTot = 0.;
363
364                t20->Coll.heat = 0.;
365                t30->Coll.heat = 0.;
366                t21->Coll.heat = 0.;
367                t31->Coll.heat = 0.;
368                t41->Coll.heat = 0.;
369
370                CoolAdd( chLabel, t20->WLAng , 0.);
371                CoolAdd( chLabel, t30->WLAng , 0.);
372                CoolAdd( chLabel, t21->WLAng , 0.);
373                CoolAdd( chLabel, t31->WLAng , 0.);
374                CoolAdd( chLabel, t41->WLAng , 0.);
375
376                /* level populations */
377                /* LIMLEVELN is the dimension of the atoms vectors */
378                ASSERT( N_SEQ_BORON <= LIMLEVELN);
379                for( i=2; i < N_SEQ_BORON; i++ )
380                {
381                        atoms.PopLevels[i] = 0.;
382                        atoms.DepLTELevels[i] = 1.;
383                }
384        }
385
386        else
387        {
388                /* atom_levelN did not evaluate PopLevels, so save pops here */
389                /* LIMLEVELN is the dimension of the atoms vectors */
390                ASSERT( N_SEQ_BORON <= LIMLEVELN);
391                for( i=0; i< N_SEQ_BORON; ++i )
392                {
393                        atoms.PopLevels[i] = pops[i];
394                        atoms.DepLTELevels[i] = depart[i];
395                }
396                /* this branch, we have a full valid solution */
397                t10->Lo->Pop = pops[0];
398                t20->Lo->Pop = pops[0];
399                t30->Lo->Pop = pops[0];
400                t21->Lo->Pop = pops[1];
401                t31->Lo->Pop = pops[1];
402                t41->Lo->Pop = pops[1];
403
404                t10->Emis->PopOpc = (pops[0] - pops[1]*t10->Lo->g/t10->Hi->g);
405                t20->Emis->PopOpc = (pops[0] - pops[2]*t20->Lo->g/t20->Hi->g);
406                t30->Emis->PopOpc = (pops[0] - pops[3]*t30->Lo->g/t30->Hi->g);
407                t21->Emis->PopOpc = (pops[1] - pops[2]*t21->Lo->g/t21->Hi->g);
408                t31->Emis->PopOpc = (pops[1] - pops[3]*t31->Lo->g/t31->Hi->g);
409                t41->Emis->PopOpc = (pops[1] - pops[4]*t41->Lo->g/t41->Hi->g);
410
411                t10->Hi->Pop = pops[1];
412                t20->Hi->Pop = pops[2];
413                t30->Hi->Pop = pops[3];
414                t21->Hi->Pop = pops[2];
415                t31->Hi->Pop = pops[3];
416                t41->Hi->Pop = pops[4];
417
418                t10->Emis->phots = t10->Emis->Aul*(t10->Emis->Pesc + t10->Emis->Pelec_esc)*pops[1];
419                t20->Emis->phots = t20->Emis->Aul*(t20->Emis->Pesc + t20->Emis->Pelec_esc)*pops[2];
420                t30->Emis->phots = t30->Emis->Aul*(t30->Emis->Pesc + t30->Emis->Pelec_esc)*pops[3];
421                t21->Emis->phots = t21->Emis->Aul*(t21->Emis->Pesc + t21->Emis->Pelec_esc)*pops[2];
422                t31->Emis->phots = t31->Emis->Aul*(t31->Emis->Pesc + t31->Emis->Pelec_esc)*pops[3];
423                t41->Emis->phots = t41->Emis->Aul*(t41->Emis->Pesc + t41->Emis->Pelec_esc)*pops[4];
424
425                t10->Emis->xIntensity = t10->Emis->phots*t10->EnergyErg;
426                t20->Emis->xIntensity = t20->Emis->phots*t20->EnergyErg;
427                t30->Emis->xIntensity = t30->Emis->phots*t30->EnergyErg;
428                t21->Emis->xIntensity = t21->Emis->phots*t21->EnergyErg;
429                t31->Emis->xIntensity = t31->Emis->phots*t31->EnergyErg;
430                t41->Emis->xIntensity = t41->Emis->phots*t41->EnergyErg;
431
432                /* ratio of collisional to total excitation */
433                t10->Emis->ColOvTot = CollRate[0][1]/(CollRate[0][1]+t10->Emis->pump);
434                t20->Emis->ColOvTot = CollRate[0][2]/(CollRate[0][2]+t20->Emis->pump);
435                t30->Emis->ColOvTot = CollRate[0][3]/(CollRate[0][3]+t30->Emis->pump);
436                t21->Emis->ColOvTot = CollRate[1][2]/(CollRate[1][2]+t21->Emis->pump);
437                t31->Emis->ColOvTot = CollRate[1][3]/(CollRate[1][3]+t31->Emis->pump);
438                t41->Emis->ColOvTot = CollRate[1][4]/(CollRate[1][4]+t41->Emis->pump);
439
440                /* derivative of cooling function */
441                thermal.dCooldT += dCoolDT;
442
443                /* two cases - collisionally excited (usual case) or
444                 * radiatively excited - in which case line can be a heat source
445                 * following are correct heat exchange, they will mix to get correct derivative
446                 * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to
447                 * keep stable solution by effectively dividing up heating and cooling,
448                 * so that negative cooling does not occur */
449
450                EnrLU = t10->Lo->Pop*CollRate[0][1]*t10->EnergyErg;
451                EnrUL = t10->Hi->Pop*CollRate[1][0]*t10->EnergyErg;
452                /* energy exchange due to this level
453                 * net cooling due to excitation minus part of de-excit */
454                t10->Coll.cool = EnrLU - EnrUL*t10->Emis->ColOvTot;
455                /* net heating is remainder */
456                t10->Coll.heat = EnrUL*(1. - t10->Emis->ColOvTot);
457                /* add to cooling */
458                CoolAdd( chLabel, t10->WLAng , t10->Coll.cool);
459                /* derivative of cooling function */
460                thermal.dCooldT += t10->Coll.cool * (t10->EnergyK * thermal.tsq1 - thermal.halfte );
461
462                EnrLU = t20->Lo->Pop*CollRate[0][2]*t20->EnergyErg;
463                EnrUL = t20->Hi->Pop*CollRate[2][0]*t20->EnergyErg;
464                t20->Coll.cool = EnrLU - EnrUL*t20->Emis->ColOvTot;
465                t20->Coll.heat = EnrUL*(1. - t20->Emis->ColOvTot);
466                /* add to cooling */
467                CoolAdd( chLabel, t20->WLAng , t20->Coll.cool);
468                thermal.dCooldT += t20->Coll.cool * (t20->EnergyK * thermal.tsq1 - thermal.halfte );
469
470                EnrLU = t30->Lo->Pop*CollRate[0][3]*t30->EnergyErg;
471                EnrUL = t30->Hi->Pop*CollRate[3][0]*t30->EnergyErg;
472                t30->Coll.cool = EnrLU - EnrUL*t30->Emis->ColOvTot;
473                t30->Coll.heat = EnrUL*(1. - t30->Emis->ColOvTot);
474                /* add to cooling */
475                CoolAdd( chLabel, t30->WLAng , t30->Coll.cool);
476                thermal.dCooldT += t30->Coll.cool * (t30->EnergyK * thermal.tsq1 - thermal.halfte );
477
478                EnrLU = t21->Lo->Pop*CollRate[1][2]*t21->EnergyErg;
479                EnrUL = t21->Hi->Pop*CollRate[2][1]*t21->EnergyErg;
480                t21->Coll.cool = EnrLU - EnrUL*t21->Emis->ColOvTot;
481                t21->Coll.heat = EnrUL*(1. - t21->Emis->ColOvTot);
482                /* add to cooling */
483                CoolAdd( chLabel, t21->WLAng , t21->Coll.cool);
484                /* use of 20 is intentional in following - that is Boltzmann factor */
485                thermal.dCooldT += t21->Coll.cool * (t20->EnergyK * thermal.tsq1 - thermal.halfte );
486
487                EnrLU = t31->Lo->Pop*CollRate[1][3]*t31->EnergyErg;
488                EnrUL = t31->Hi->Pop*CollRate[3][1]*t31->EnergyErg;
489                t31->Coll.cool = EnrLU - EnrUL*t31->Emis->ColOvTot;
490                t31->Coll.heat = EnrUL*(1. - t31->Emis->ColOvTot);
491                /* add to cooling */
492                CoolAdd( chLabel, t31->WLAng , t31->Coll.cool);
493                /* use of 30 is intentional in following - that is Boltzmann factor */
494                thermal.dCooldT += t31->Coll.cool * (t30->EnergyK * thermal.tsq1 - thermal.halfte );
495
496                EnrLU = t41->Lo->Pop*CollRate[1][4]*t41->EnergyErg;
497                EnrUL = t41->Hi->Pop*CollRate[4][1]*t41->EnergyErg;
498                t41->Coll.cool = EnrLU - EnrUL*t41->Emis->ColOvTot;
499                t41->Coll.heat = EnrUL*(1. - t41->Emis->ColOvTot);
500                /* add to cooling */
501                CoolAdd( chLabel, t41->WLAng , t41->Coll.cool);
502                /* use of 41 is intentional in following - that is Boltzmann factor (no 40 here) */
503                thermal.dCooldT += t41->Coll.cool * (t41->EnergyK * thermal.tsq1 - thermal.halfte );
504        }
505
506        return;
507}
Note: See TracBrowser for help on using the browser.