cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_co_atom.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 /*CO_PopsEmisCool evaluate rotation levels populations, emission, and cooling */
4 /*cdCO_colden return column density in H2, negative -1 if cannot find state,
5  * header is cddrive */
6 #include "cddefines.h"
7 #include "taulines.h"
8 #include "dense.h"
9 #include "thermal.h"
10 #include "hmi.h"
11 #include "radius.h"
12 #include "atoms.h"
13 #include "phycon.h"
14 #include "rt.h"
15 #include "lines_service.h"
16 #include "cddrive.h"
17 #include "mole.h"
18 /*lint -e778 constant expression evaluatess to 0 in operation '-' */
19 /*=================================================================*/
20 
21 /* will be used to save CO column densities */
22 static double *col12 , *col13;
23 
24 /* evaluate CO rotation cooling */
26  transition ** Rotate ,
27  long int nRotate ,
28  realnum abundan,
29  const char * chLabel ,
30  realnum * Cooling ,
31  realnum * dCoolingDT )
32 {
33 
34  /* will need to MALLOC space for these but only on first call */
35  static double
36  **AulEscp ,
37  **col_str ,
38  **AulDest,
39  /* AulPump[low][high] is rate (s^-1) from lower to upper level */
40  **AulPump,
41  **CollRate,
42  *pops,
43  *create,
44  *destroy,
45  *depart,
46  /* statistical weight */
47  *stat ,
48  /* excitation energies in kelvin */
49  *excit;
50 
51  /*static long int **ipdest;*/
52 
53  static bool lgFirst=true;
54  long int i,
55  j,
56  ilo ,
57  ihi;
58  static long int nUsed;
59  bool lgDeBug;
60  int nNegPop;
61  bool lgZeroPop;
62  double rot_cooling , dCoolDT;
63  static long int ndimMalloced = 0;
64 
65  DEBUG_ENTRY( "CO_PopsEmisCool()" );
66 
67  if( lgFirst )
68  {
69  /* will never do this again */
70  lgFirst = false;
71  /* remember how much space we malloced in case ever called with more needed */
72  ndimMalloced = nRotate;
73  /* allocate the 1D arrays*/
74  excit = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
75  stat = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
76  pops = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
77  create = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
78  destroy = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
79  depart = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
80  /* create space for the 2D arrays */
81  AulPump = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
82  CollRate = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
83  AulDest = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
84  AulEscp = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
85  col_str = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
86  for( i=0; i<(nRotate+1); ++i )
87  {
88  AulPump[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
89  CollRate[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
90  AulDest[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
91  AulEscp[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
92  col_str[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
93  }
94  }
95  /* this is test for call with too many rotation levels to handle - logic needs
96  * for largest rotor to be called first */
97  if( nRotate > ndimMalloced )
98  {
99  fprintf(ioQQQ," CO_PopsEmisCool has been called with the number of rotor levels greater than space allocated.\n");
100  cdEXIT(EXIT_FAILURE);
101  }
102 
103  /* all elements are used, and must be set to zero if zero */
104  for( i=0; i < (nRotate+1); i++ )
105  {
106  create[i] = 0.;
107  destroy[i] = 0.;
108  for( j=0; j < (nRotate+1); j++ )
109  {
110  col_str[j][i] = 0.;
111  AulEscp[j][i] = 0.;
112  AulDest[j][i] = 0.;
113  AulPump[j][i] = 0.;
114  }
115  }
116 
117  /* the statistical weights of the levels */
118  for( j=0; j < nRotate; j++ )
119  {
120  /* statistical weights for each level */
121  stat[j] = (*Rotate)[j].Lo->g;
122  }
123  /* this is the highest level, which is one more than the highest line */
124  stat[nRotate] = (*Rotate)[nRotate-1].Hi->g;
125 
126  /* set up the excitation potentials of each level relative to ground -
127  * the struc saves the energy of the line only */
128  excit[0] = 0.;
129  for( j=1; j < nRotate; j++ )
130  {
131  /* excitation energy of each level relative to ground, in K */
132  excit[j] = excit[j-1] + (*Rotate)[j-1].EnergyK;
133  }
134 
135  /* this is the highest level, which is one more than the highest line */
136  excit[nRotate] = excit[nRotate-1] + (*Rotate)[nRotate-1].EnergyK;
137 
138  nUsed = nRotate;
139 
140  /* this determines the largest molecule that can be inverted at this
141  * temperature. Need Boltzmann factors to be positive for all levels
142  * nUsed is the index of the highest level, so the number of levels (passed
143  * to the solver) is nUsed+1
144  * excit[nUsed]/phycon.te cannot be much larger than 20, or matrix inversion will fail */
145  /* >>chng 03 sep 18, allow to go to 1, a two level atom */
146  /*while ( (excit[nUsed] > phycon.te*20.) && (nUsed > 1) )*/
147  /* >>chng 03 oct 03, chng 20 to 25 */
148  /*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/
149  /* >>chng 03 oct 03, keep at least 5 levels */
150  /*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/
151  while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 5) )
152  {
153  --nUsed;
154  }
155 
156  for( j=0; j < nRotate; j++ )
157  {
158  /*data[j][j+1] = (*Rotate)[j].Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc);*/
159  AulEscp[j+1][j] = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc);
160  AulDest[j+1][j] = (*Rotate)[j].Emis->Aul*(*Rotate)[j].Emis->Pdest;
161  /* next 2 not not used by levelN since flag passed saying to use CollRate instead */
162  (*Rotate)[j].Coll.cs = 1.;
163  /*data[j+1][j] = (*Rotate)[j].Coll.cs*hmi.Hmolec[ipMH2g]/dense.eden;*/
164  col_str[j+1][j] = (*Rotate)[j].Coll.cs*hmi.H2_total/dense.eden;
165  /* photon pumping rate */
166  AulPump[j][j+1] = (*Rotate)[j].Emis->pump;
167  /* the continuum indices are on the f, not c, scale, and will be passed to
168  * atom_levelN, which works on f, not c, scale */
169  /*ipdest[j][j+1] = (*Rotate)[j].ipCont;*/
170  }
171 
172  /* next loop for collisional transitions that have delta J >1 */
173  for( ilo=0; ilo < nRotate-1; ilo++ )
174  {
175  /* need to have upper limit to to nRotate because
176  * number ov levels is 1 gt number of lines*/
177  for( ihi=ilo+2; ihi <= nRotate; ihi++ )
178  {
179  /* this is not used by levelN since flag passed saying to use CollRate instead */
180  /*data[ihi][ilo] = 1. *hmi.Hmolec[ipMH2g]/dense.eden;*/
181  /*data[ihi][ilo] = 1. *hmi.H2_total/dense.eden;*/
182  col_str[ihi][ilo] = 1. *hmi.H2_total/dense.eden;
183  /* these are escape and dest rates, which are zero for a rigid rotator */
184  /*data[ilo][ihi] = 0;*/
185  AulEscp[ihi][ilo] = 0;
186  AulDest[ihi][ilo] = 0;
187  }
188  }
189 
190  /* now evaluate the H2 collision rates */
191  /* recall one more level than lines */
192  for( ilo=0; ilo < nRotate; ilo++ )
193  {
194  /* need to have upper limit to to nRotate because
195  * number of levels is 1 gt number of lines*/
196  for( ihi=ilo+1; ihi <= MIN2(nRotate,ilo+5); ihi++ )
197  {
198  /* >>refer CO-H2 collision de Jong, T., Chu, S-I., & Dalgarno, A. 1975, ApJ, 199, 69 */
199  double a[5]={1.66, 2.80, 1.19, 1.00, 1.30};
200  double b[5]={1.67, 1.47, 1.85, 1.55, 2.24};
201  /* >>refer CO-He collision McKee, C.F., Storey, J.W.V., Watson, D.M., & Green, S., 1982, ApJ, 259, 647 */
202  /* this reference gives He-CO collisions,
203  * their table 1 says He collisions are 1.37 slower than H2 collisions*/
204  /* >>chng 03 sep 19, from ground H2 to total H2 */
205  double collid = hmi.H2_total + dense.xIonDense[ipHELIUM][0]/1.37;
206  /* de Jong et al. explicitly consider temperatures as low as 20K,
207  * don't extrapolate significantly lower than this */
208  double TeUsed = MAX2(10., phycon.te );
209 
210  /* first do deexcitation rate, equation 17 of deJong et al. */
211  CollRate[ihi][ilo] = a[ihi-ilo-1]*1.e-10*(*Rotate)[ilo].Lo->g/(*Rotate)[ilo].Hi->g*
212  (1.+(excit[ihi]-excit[ilo])/TeUsed) *
213  sexp( b[ihi-ilo-1]*sqrt((excit[ihi]-excit[ilo])/TeUsed) )*collid;
214  /* this is mainly so that rest of code gets valid cs in case crit dense printed */
215  if( ihi == ilo+1 )
216  {
217  /* save downward collision rate */
218  (*Rotate)[ilo].Coll.ColUL = (realnum)CollRate[ihi][ilo];
219  /* convert into fake electron collision strength */
220  LineConvRate2CS( &(*Rotate)[ilo] , (*Rotate)[ilo].Coll.ColUL);
221  }
222 
223  /* now get excitation rate */
224  CollRate[ilo][ihi] = CollRate[ihi][ilo]*
225  sexp( (excit[ihi]-excit[ilo])/phycon.te )*
226  (*Rotate)[ilo].Hi->g / (*Rotate)[ilo].Lo->g;
227 
228  /* debug print statement */
229  /*fprintf(ioQQQ," %li %li %.2e %.2e \n",ilo, ihi,
230  CollRate[ihi][ilo]/hmi.H2_total , CollRate[ilo][ihi]/hmi.H2_total );*/
231  }
232  /* finish off with zeros */
233  for( ihi=ilo+6; ihi <= nRotate; ihi++ )
234  {
235  CollRate[ihi][ilo] = 0.;
236  CollRate[ilo][ihi] = 0.;
237  }
238  }
239 
240  lgDeBug = false;
241 
242  atom_levelN(
243  /* number of levels is number of lines plus one */
244  /* set nUsed so that CO is evaluated even at very low temperatures */
245  nUsed+1, /*nRotate+1,*/
246  abundan,
247  stat,
248  excit,
249  'K',
250  pops,
251  depart,
252  &AulEscp,
253  &col_str,
254  &AulDest,
255  &AulPump,
256  &CollRate,
257  create,
258  destroy,
259  /* say that we have evaluated the collision rates already */
260  true,
261  /*&ipdest,*/
262  &rot_cooling,
263  &dCoolDT,
264  chLabel,
265  /* nNegPop positive if negative pops occured, negative if too cold */
266  &nNegPop,
267  &lgZeroPop,
268  lgDeBug );/* option to print stuff - set to true for debug printout */
269 
270  /* put cooling into place where we can use it later */
271  *Cooling = (realnum)rot_cooling;
272 
273  /* >>chng 03 sep 18, CO rot cooling temp deriv is too small, and temp
274  * changes too big - incr deriv by 5x
275  *dCoolingDT = (realnum)(dCoolDT);*/
276  *dCoolingDT = (realnum)(rot_cooling/phycon.te);
277 
278 # if 0
279  if( lgMainCO )
280  fprintf(ioQQQ,"COcool\t%i\t%.2f\t%e\t%e\t%e\t%e\t%e\n",
281  nUsed,
282  fnzone,
283  phycon.te,
284  *Cooling/abundan,
285  *dCoolingDT ,
286  *Cooling,
287  thermal.htot );
288 # endif
289 
290  /* zero out higher populations for case where full CO levels not done */
291  for( i=nUsed+1; i<=nRotate; ++i )
292  {
293  pops[i] = 0.;
294  depart[i] = 0.;
295  }
296 
297  /* can only define first LIMLEVELN elements, the vector's length */
298  for( i=0; i< MIN2(LIMLEVELN,nRotate); ++i )
299  {
300  atoms.PopLevels[i] = pops[i];
301  atoms.DepLTELevels[i] = depart[i];
302  }
303 
304  if( nNegPop > 0 )
305  {
306  fprintf(ioQQQ,"CO_PopsEmisCool called atom_levelN which returned negative populations.\n");
307  }
308 
309  /* now establish information that is passed out to rest of code's infrastructure */
310  for( j=0; j<nRotate; ++j )
311  {
312  double EnrLU, EnrUL;
313  /* lower upper populations, stim emission correct population */
314  (*Rotate)[j].Lo->Pop = pops[j];
315  (*Rotate)[j].Hi->Pop = pops[j+1];
316  (*Rotate)[j].Emis->PopOpc = (pops[j] - pops[j+1]*(*Rotate)[j].Lo->g/(*Rotate)[j].Hi->g);
317 
318  /* number of photons in the line */
319  (*Rotate)[j].Emis->phots = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc)*pops[j+1];
320 
321 # if 0
322  /* >>chng 03 oct 04, moved to RT_OTS */
323  /* local Emis.ots rates, added to line ots array */
324  (*Rotate)[j].Emis->ots = (*Rotate)[j].Aul*(*Rotate)[j].Emis->Pdest*(realnum)(*Rotate)[j].Hi->Pop;
325  RT_OTS_AddLine( (*Rotate)[j].Emis->ots , (*Rotate)[j].ipCont);
326 # endif
327 
328  /* the intensity in the line */
329  (*Rotate)[j].Emis->xIntensity = (*Rotate)[j].Emis->phots*(*Rotate)[j].EnergyErg;
330 
331  /* ratio of collisional to total excitation */
332  (*Rotate)[j].Emis->ColOvTot = (realnum)(CollRate[j][j+1]/(CollRate[j][j+1]+(*Rotate)[j].Emis->pump) );
333 
334  /* two cases - collisionally excited (usual case) or
335  * radiatively excited - in which case line can be a heat source
336  * following are correct heat exchange, if will mix to get correct deriv */
337  EnrLU = (*Rotate)[j].Lo->Pop*CollRate[j][j+1]*(*Rotate)[j].EnergyErg;
338  EnrUL = (*Rotate)[j].Hi->Pop*CollRate[j+1][j]*(*Rotate)[j].EnergyErg;
339  /* energy exchange due to this level
340  * net cooling due to excit minus part of de-excit */
341  (*Rotate)[j].Coll.cool = EnrLU - EnrUL*(*Rotate)[j].Emis->ColOvTot;
342  /* net heating is remainder */
343  (*Rotate)[j].Coll.heat = EnrUL*(1.f - (*Rotate)[j].Emis->ColOvTot);
344  /* do not add to cooling, since done with evaluated cooling from atom_levelN */
345  /*CoolAdd( chLabel, (long)t10->WLAng , t10->cool);*/
346  }
347 
348  /* generate flag if co cooling important and highest level is fainter
349  * than second highest level */
350  if( rot_cooling / thermal.ctot > 0.1 &&
351  (*Rotate)[nUsed-1].Emis->xIntensity > (*Rotate)[nUsed-2].Emis->xIntensity &&
352  /*>>chng 03 oct 03, add check whether molecule has been trimed down
353  * due to small Boltzmann factor */
354  /* >>chng 03 oct 20, following had; after ), so CoolCaped always true */
355  nUsed == nRotate )
356  {
357  co.lgCOCoolCaped = true;
358  }
359 
360  return;
361 }
362 
363 /*CO_Colden maintain H2 column densities within X */
364 void CO_Colden( const char *chLabel )
365 {
366  long int iRot;
367 
368  DEBUG_ENTRY( "CO_Colden()" );
369 
370  if( strcmp(chLabel,"ZERO") == 0 )
371  {
372  static bool lgFIRST=true;
373  if( lgFIRST )
374  {
375  lgFIRST = false;
376  col12 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) );
377  col13 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) );
378  }
379  /* the column density (cm-2) of ortho and para H2 */
380  /* zero out formation rates and column densites */
381  for( iRot=0; iRot<=nCORotate; ++iRot )
382  {
383  /* zero it out */
384  col12[iRot] = 0.;
385  col13[iRot] = 0.;
386  }
387  }
388  else if( strcmp(chLabel,"ADD ") == 0 )
389  {
390  /* add together column densities */
391  for( iRot=0; iRot<nCORotate; ++iRot )
392  {
393  /* zero it out */
394  col12[iRot] += C12O16Rotate[iRot].Lo->Pop*radius.drad_x_fillfac;
395  col13[iRot] += C13O16Rotate[iRot].Lo->Pop*radius.drad_x_fillfac;
396  }
397  }
398 
399  /* we will not print column densities so skip that - if not print then we have a problem */
400  else if( strcmp(chLabel,"PRIN") != 0 )
401  {
402  fprintf( ioQQQ, " CO_Colden does not understand the label %s\n",
403  chLabel );
404  cdEXIT(EXIT_FAILURE);
405  }
406 }
407 
408 /*cdCO_colden return column density in H2, negative -1 if cannot find state,
409  * header is cddrive */
410 double cdCO_colden( long isotope , long iRot )
411 {
412 
413  /* make sure incoming parameters are ok */
414  if( isotope !=12 && isotope != 13 )
415  {
416  fprintf(ioQQQ," cdCO_colden can only do 12CO and 13CO\n");
417  return -1.;
418  }
419  if( iRot < 0 || iRot > nCORotate )
420  {
421  fprintf(ioQQQ," cdCO_colden - rotation quantum number must be 0 or greater, and less than %li\n",
422  nCORotate);
423  return -1.;
424  }
425 
426  /* the incoming parameters are fully qualified - return the column density */
427  if( isotope == 12 )
428  {
429  return col12[iRot];
430  }
431  else
432  {
433  return col13[iRot];
434  }
435 }
436 
437 /*CO_OTS - add CO lines to ots fields */
438 void CO_OTS( void )
439 {
440 
441  long int j;
442 
443  DEBUG_ENTRY( "CO_OTS()" );
444 
445  /* add all CO lines */
446  for( j=0; j<nCORotate; ++j )
447  {
449  C12O16Rotate[j].Hi->Pop;
450  RT_OTS_AddLine( C12O16Rotate[j].Emis->ots , C12O16Rotate[j].ipCont);
451 
453  C13O16Rotate[j].Hi->Pop;
454  RT_OTS_AddLine( C13O16Rotate[j].Emis->ots , C13O16Rotate[j].ipCont);
455  }
456 
457  return;
458 }

Generated for cloudy by doxygen 1.8.3.1