cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cool_nitr.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 /*CoolNitr evaluate total cooling due to nitrogen */
4 #include "cddefines.h"
5 #include "taulines.h"
6 #include "dense.h"
7 #include "phycon.h"
8 #include "ligbar.h"
9 #include "thermal.h"
10 #include "lines_service.h"
11 #include "embesq.h"
12 #include "atoms.h"
13 #include "cooling.h"
14 #include "nitro.h"
15 /* xNI_coll_stren: Author: Terry Yun 2006 */
16 /* This program interpolate the collision strength of NI
17  * We are only interested in first 5 levels, therefor there are 10 transitions.
18  * Used the Excel to find the polynomial fit, eqn: CS = a*t+b*t^1.5+c*t^0.5
19  * The range of temp is 0 - 50,000 K */
20 /* #define PRINT */
21 static const int N1_SIZE = 10;
22 
23 static double aNI[N1_SIZE] = {2.755e-5,4.123e-5,7.536e-6,1.486e-5,4.516e-5,-2.935e-6,4.000e-6,3.751e-6,-2.176e-6,1.024e-5};
24 static double bNI[N1_SIZE] = {-8.150e-8,-1.220e-7,-2.226e-8,-4.390e-8,-1.130e-7,8.000e-9,-1.1447e-8,-1.061e-8,5.610e-9,-3.227e-8};
25 static double cNI[N1_SIZE] = {2.140e-4,3.272e-4,-4.944e-6,3.473e-6,-8.772e-4,1.654e-3,1.675e-3,1.123e-3,3.867e-3,3.376e-4};
26 
27 static double NI[6][6][3]; /* coefficients a b c into array */
28 
29 /*xNI_coll_stren interpolate the collision strength of NI */
30 STATIC double xNI_coll_stren(int init, int final)
31 {
32  /* incorporating the N0 collision strengths give in
33  *>>refer N1 cs Tayal, S.S., 2006, ApJS, 163, 207 */
34 
35  int i, j;
36  double CS; /* collision strength */
37  int index = 0;
38  double temp_max = 50e3;
39  double temp_min = 0;
40  double temp = 0;
41 
42  /* assigning array initially to zero */
43  for(i=0; i<6; i++)
44  {
45  for(j=0; j<6; j++)
46  {
47  NI[i][j][0] = 0; /* coeff a */
48  NI[i][j][1] = 0; /* coeff b */
49  NI[i][j][2] = 0; /* coeff c */
50  }
51  }
52 
53  /* reading coeffs into 3D array */
54  /* The index is in physics scale */
55  for(i=1; i<6; i++)
56  {
57  for(j=i+1; j<6; j++)
58  {
59  NI[i][j][0] = aNI[index];
60  NI[i][j][1] = bNI[index];
61  NI[i][j][2] = cNI[index];
62  index++;
63  }
64  }
65 
66 # ifdef PRINT
67  /* physical index goes from 1 to 5 */
68  for(i=1; i<6; i++)
69  {
70  for(j=i+1; j<6; j++)
71  {
72  printf("NI %i%i a:%e ", i, j, NI[i][j][0]);
73  printf("%i%i b:%e\n", i, j, NI[i][j][1]);
74  /* printf("%i%i c:%lf\n", i, j, NI[i][j][2]); */
75  }
76  }
77 # endif
78 
79 
80  /* Invalid entries returns '-1':the initial indices are smaller than the final indices */
81  if(init >= final)
82  {
83  CS = -1;
84  }
85 
86  /* Invalid returns '-1': the indices are greater than 5 or smaller than 0 */
87  else if(init < 1 || init > 5 || final < 1 || final > 5)
88  {
89  CS = -1;
90  }
91 
92  else
93  {
94  temp = MAX2(temp, temp_min); /* return lager temp */
95  temp = MIN2(temp, temp_max); /* return smaller temp */
96  /* eqn: CS = a*t+b*t^1.5+c*t^0.5 */
97  CS = NI[init][final][0]*phycon.te + NI[init][final][1]*phycon.te32 + NI[init][final][2]*phycon.sqrte;
98  }
99 
100  return CS;
101 }
102 
103 
104 void CoolNitr(void)
105 {
106  realnum p2;
107  double a21,
108  a31,
109  a32,
110  cs,
111  cs2s2p,
112  cs2s3p,
113  cs21,
114  cs31,
115  cs32,
116  p3,
117  pump_rate;
118  long int i;
119 
120  double a12,
121  a13,
122  a14,
123  a15,
124  a23,
125  a24,
126  a25,
127  a34,
128  a35,
129  a45,
130  cs12,
131  cs13,
132  cs14,
133  cs15,
134  cs23,
135  cs24,
136  cs25,
137  cs34,
138  cs35,
139  cs45;
140 
141  double pop[5];
142  static double gAr4[5]={4.,6.,4.,2.,4.};
143  static double exAr4[4]={19224.464,8.713,9605.74,0.386};
144 
145  static long int *ipN1Pump=NULL,
146  nN1Pump=0 , lgFirst_N1=true;
147 
148  static bool lgFirst_N3=true;
149  static long int *ipN3Pump=NULL,
150  nN3Pump=0;
151 
152 /* #define PRINT */
153 
154 # ifdef PRINT
155  FILE *ofp;
156  double pop1, pop2,pop3,pop4,pop5;
157  ofp = fopen("c:\\projects\\cloudy\\trunk\\mydata\\pop.out", "w");
158  if(ofp == NULL)
159  {
160  printf("Can't open the file\n");
161  cdEXIT( EXIT_FAILURE );
162  }
163 # endif
164 
165  DEBUG_ENTRY( "CoolNitr()" );
166 
167  /* >>chng 00 dec 11, add made-up collision strengths
168  * These have made-up collision strengths */
169 # if 0
170  PutCS(1.,&TauLines[ipT671]);
171  atom_level2(&TauLines[ipT671]);
172 
173  PutCS(1.,&TauLines[ipT315]);
174  atom_level2(&TauLines[ipT315]);
175 
176  PutCS(1.,&TauLines[ipT333]);
177  atom_level2(&TauLines[ipT333]);
178 
179  PutCS(1.,&TauLines[ipT324]);
180  atom_level2(&TauLines[ipT324]);
181 
182  PutCS(1.,&TauLines[ipT374g]);
183  atom_level2(&TauLines[ipT374g]);
184 
185  PutCS(1.,&TauLines[ipT374x]);
186  atom_level2(&TauLines[ipT374x]);
187 # endif
188 
189  /* find pump rate for [N I] 5199, do one time initialization if first
190  * call and level 2 lines are on - NB lgFirst is used in two places
191  * in this routine - here we use it, and the last place where it is
192  * used it will be set false */
193  if( lgFirst_N1 && nWindLine && dense.lgElmtOn[ipNITROGEN] )
194  {
195  lgFirst_N1 = false;
196  nN1Pump = 0;
197  for( i=0; i<nWindLine; ++i )
198  {
199  /* don't test on nelem==ipNITROGEN since lines on physics, not C, scale */
200  if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==1 )
201  {
202  ++nN1Pump;
203  }
204  }
205  if( nN1Pump<0 )
206  TotalInsanity();
207  else if( nN1Pump > 0 )
208  /* create the space - can't malloc 0 bytes */
209  ipN1Pump = (long *)MALLOC((unsigned)(nN1Pump)*sizeof(long) );
210  nN1Pump = 0;
211  for( i=0; i<nWindLine; ++i )
212  {
213  /* don't test on nelem==ipNITROGEN since lines on physics, not C, scale */
214  if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==1 )
215  {
216 # if 0
217  DumpLine( &TauLine2[i] );
218 # endif
219  ipN1Pump[nN1Pump] = i;
220  ++nN1Pump;
221  }
222  }
223  }
224  else
225  /* level 2 lines are not enabled */
226  nN1Pump = 0;
227 
228  /* now sum pump rates - will be used in prt statement together
229  * with quench rate derive here to estimate continuum pumping of
230  [N I] 5199 */
231  nitro.pump_rate_N1 = 0.;
232  for( i=0; i<nN1Pump; ++i )
233  {
234  nitro.pump_rate_N1 += TauLine2[ipN1Pump[i]].Emis->pump;
235 # if 0
236  fprintf(ioQQQ,"DEBUG N %li %.3e %.3e\n",
237  i,
238  TauLine2[ipN1Pump[i]].WLAng , TauLine2[ipN1Pump[i]].pump );
239 # endif
240  }
241 
242  /* >>chng 01 sep 02, go back to old values */
243  /* nitrogen I 5200, cs and A from
244  * collision strengths from
245  * >>referold n1 cs Dopita, M.A., Mason, D.J., Robb, W.D. 1976, ApJ, 207, 102*/
246  cs21 = 1.32e-4*phycon.te/(phycon.te10*phycon.te01);
247  cs31 = 3.60e-5*phycon.te/phycon.te10*phycon.te02;
248  cs32 = 5.23e-4*phycon.te70*phycon.te10/phycon.te02;
249 
250  /* N.B. following cs also appears in nitrogen ionization balance */
251  /* >>chng revise collision strengths to
252  * >>refer n1 cs Tayal, S.S., 2000, ADNDT, 76, 191
253  cs21 = 5.25e-13*phycon.tesqrd*phycon.te70*phycon.te03;
254  cs31 = 3.79e-5*phycon.te70*phycon.te03;
255  cs32 = 3.14e-13*phycon.tesqrd*phycon.te*phycon.te20*phycon.te07/phycon.te003;*/
256 
257  /* following does not have photoionization from excited state
258  * since already counted in photo balance equation */
261  /* trans prob from
262  * >>refer n1 as Butler, K., & Zeippen, C.J. 1984, A&A, 141, 274
263  */
264  a21 = 1.28e-5;
265  a31 = 5.31e-3;
266  a32 = 6.81e-2;
267  /* >>chng 01 sep 02, include atoms.d5200r term here too */
268 
271  /* trans prob from
272  * >>refer n1 as Butler, K., & Zeippen, C.J. 1984, A&A, 141, 274
273  */
274  cs12 = xNI_coll_stren(1, 2);
275  a12 = 7.170e-6;
276 
277  cs13 = xNI_coll_stren(1, 3);
278  a13 = 3.994e-5;
279 
280  cs14 = xNI_coll_stren(1, 4);
281  a14 = 5.397e-3;
282 
283  cs15 = xNI_coll_stren(1,5);
284  a15 = 1.329e-2;
285 
286  /* this is fraction of 2D excitations that emit a photon in the 5200 multiplet
287  * this is used for collisional quenching in prt lines routine */
288  nitro.quench_5200 = (a12+a13) / ((a12+a13) + (cs12 + cs13) * dense.cdsqte );
289 
290  /* collision strengths and transition probabilities for the lowest 5
291  * levels of N^0 */
292  cs23 = xNI_coll_stren(2, 3);
293  a23 = 2.478e-8;
294 
295  cs24 = xNI_coll_stren(2, 4);
296  a24 = 3.135e-2;
297 
298  cs25 = xNI_coll_stren(2, 5);
299  a25 = 5.693e-2;
300 
301  cs34 = xNI_coll_stren(3, 4);
302  a34 = 4.92e-2;
303 
304  cs35 = xNI_coll_stren(3, 5);
305  a35 = 2.709e-2;
306 
307  cs45 = xNI_coll_stren(4, 5);
308  a45 = 1e-15;
309 
310  atom_pop5(gAr4,exAr4,cs12,cs13,cs14,cs15,cs23,cs24,cs25,cs34,cs35,
311  cs45,a12,a13,a14,a15,a23,a24,a25,a34,a35,a45,pop,dense.xIonDense[ipNITROGEN][0]);
312 
313 # ifdef PRINT
314  pop1 = pop[0]/dense.xIonDense[ipNITROGEN][0];
315  pop2 = pop[1]/dense.xIonDense[ipNITROGEN][0];
316  pop3 = pop[2]/dense.xIonDense[ipNITROGEN][0];
317  pop4 = pop[3]/dense.xIonDense[ipNITROGEN][0];
318  pop5 = pop[4]/dense.xIonDense[ipNITROGEN][0];
319  fprintf(ofp, "%lf %lf %lf %lf %lf\n", pop1, pop2, pop3, pop4, pop5);
320 
321 # endif
322 
323  nitro.xN5200 = pop[1]*a12*3.818e-12;
324  nitro.xN5198 = pop[2]*a13*3.821e-12;
325  nitro.xN3467 = pop[3]*a14*5.729e-12;
326  nitro.xN3466 = pop[4]*a15*5.729e-12;
327  nitro.xN10398 = pop[3]*a24*1.910e-12;
328  nitro.xN10397 = pop[4]*a25*1.910e-12;
329  nitro.xN10408 = pop[3]*a34*1.908e-12;
330  nitro.xN10407 = pop[4]*a35*1.908e-12;
331 
332  p3 = atom_pop3(4.,10.,6.,cs21,cs31,cs32,a21,a31,a32,2.769e4,1.38e4,
333  &p2,dense.xIonDense[ipNITROGEN][0],atoms.d5200r,0.,0.);
334  /* save population */
336 
337  nitro.c5200 = p2*a21*3.83e-12;
338  CoolAdd("N 1",5200,nitro.c5200);
340 
341  nitro.c10400 = p3*a32*1.91e-12;
342  CoolAdd("N 1",10400,nitro.c10400);
343 
344  nitro.c3466 = p3*a31*5.74e-12;
345  CoolAdd("N 1",3466,nitro.c3466);
346  thermal.dCooldT += (nitro.c10400 + nitro.c3466)*(4.15e4*
348 
349  /* N I 1200, cs from trans probability */
350  PutCS(4.1,&TauLines[ipT1200]);
351  atom_level2(&TauLines[ipT1200]);
352 
353  /* N II 1084, cs from trans prob */
354  PutCS(5.5,&TauLines[ipT1085]);
355  atom_level2(&TauLines[ipT1085]);
356 
357  /* [N II] coll data from
358  * >>referold n2 cs Stafford, R.P., Bell, K.L, Hibbert, A. & Wijesundera, W.P.,
359  * >>rereroldcon 1994, MNRAS 268, 816,
360  * at 10,000K (v weak T dep)
361  * >>chng 00 dec 11, to Lennon & Burke
362  * >>refer n2 cs Lennon, D.J., & Burke, V.M., 1994, A&AS, 103, 273-277
363  * transit prob from
364  * >>refer n2 as Nussbaumer, H., & Rusca, C. 1979, A&A, 72, 129
365  *
366  * >>chng 00 dec 11, to Lennon & Burke cs */
367  a21 = 3.65e-3;
368  a31 = 0.0316;
369  a32 = 1.17;
370  /* >>chng 02 may 02, put in option to switch between Lennon Burke and Stafford et al. ,
371  * default state of this var is true */
375 # define USE_LENNON_BURKE 1
376 # if( USE_LENNON_BURKE )
377  cs21 = 2.64;
378  cs31 = 0.293;
379  cs32 = 0.834;
380 # else
381  /* Bob Rubin could not find the Stafford et al. numbers - it flag ever
382  * set false will not compile */
383  cs21 = 3.02;
384  cs31 = 0.372;
385  cs32 = 0.505;
386 # endif
387 
388  /* POP3( G1,G2,G3, O12, O13, O23, A21, A31, A32,*/
389  p3 = atom_pop3( 9.,5.,1., cs21, cs31, cs32, a21, a31, a32 ,
390  /* E12,E23,P2,ABUND,GAM2) */
391  21955.,24982.,&p2,dense.xIonDense[ipNITROGEN][1],0.,0.,0.);
392 
393  nitro.c5755 = p3*a32*3.46e-12;
394 
395  nitro.c6584 = p2*a21*3.03e-12;
397  CoolAdd("N 2",5755,nitro.c5755);
398  CoolAdd("N 2",6584,nitro.c6584);
399 
400  /* nitro.xN2_A3_tot is fraction of excitations that produce a photon
401  * and represents the correction for collisiona deexcitation */
402  nitro.xN2_A3_tot = (a31+a32) /(a31+a32 + (cs31+cs32)/1.*dense.cdsqte );
403  ASSERT( nitro.xN2_A3_tot <= 1. );
404 
405  /* N II fine structure lines, */
406  /* >>chng 00 dec 11, to Lennon & Burke CS */
407 # if( USE_LENNON_BURKE )
408  cs21 = 0.408;
409  cs32 = 1.12;
410  cs31 = 0.272;
411 # else
412  /* Bob Rubin could not find the Stafford et al. numbers - it flag ever
413  * set false will not compile */
414  cs21 = 0.429;
415  cs32 = 1.13;
416  cs31 = 0.265;
417 # endif
418 
419  PutCS(cs21,&TauLines[ipT205]);
420  PutCS(cs32,&TauLines[ipT122]);
421  PutCS(cs31,&TauDummy);
422  atom_level3(&TauLines[ipT205],&TauLines[ipT122],&TauDummy);
423 
424  /* N II 2140, data
425  * >>refer n2 cs Stafford, R.P., Bell, K.L, Hibbert, A. & Wijesundera, W.P. 1994,
426  * >>rerercon MNRAS, 268, 816
427  * >>refer n2 cs Lennon, D.J., & Burke, V.M., 1994, A&AS, 103, 273-277
428  * A from
429  * >>refer n2 as Brage, T., Hibbert, A., Leckrone, D.S. 1997, ApJ, 478, 423 */
430 # if( USE_LENNON_BURKE )
431  cs21 = 1.19;
432 # else
433  /* Bob Rubin could not find the Stafford et al. numbers - it flag ever
434  * set false will not compile */
435  cs21 = 1.15;
436 # endif
437 
438  PutCS(cs21,&TauLines[ipT2140]);
439  atom_level2(&TauLines[ipT2140]);
440 
441  /* N III 989.8, cs from
442  * >>refer N3 cs Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 */
443  PutCS(7.12,&TauLines[ipT990]);
444  atom_level2(&TauLines[ipT990]);
445 
446  /* 57 micron N III, A=
447  * >>refer n3 as Froese Fischer, C. 1983, J.Phys. B, 16, 157
448  * collision strength from
449  * >>refer n3 cs Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 */
450  cs = MIN2(1.90,0.2291*phycon.te10*phycon.te10);
451  PutCS(cs,&TauLines[ipT57]);
452 
453 
454  /* one time initialization if first call, and level 2 lines are on */
455  if( lgFirst_N3 && nWindLine && dense.lgElmtOn[ipNITROGEN] )
456  {
457  lgFirst_N3 = false;
458  nN3Pump = 0;
459  for( i=0; i<nWindLine; ++i )
460  {
461  /* don't test on nelem==ipIRON since lines on physics, not C, scale */
462  if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==3 )
463  {
464  ++nN3Pump;
465  }
466  }
467  if( nN3Pump<0 )
468  TotalInsanity();
469  else if( nN3Pump > 0 )
470  /* create the space - can't malloc 0 bytes */
471  ipN3Pump = (long *)MALLOC((unsigned)(nN3Pump)*sizeof(long) );
472  nN3Pump = 0;
473  for( i=0; i<nWindLine; ++i )
474  {
475  /* don't test on nelem==ipIRON since lines on physics, not C, scale */
476  if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==3 )
477  {
478 # if 0
479  DumpLine( &TauLine2[i] );
480 # endif
481  ipN3Pump[nN3Pump] = i;
482  ++nN3Pump;
483  }
484  }
485  }
486  else
487  /* level 2 lines are not enabled */
488  nN3Pump = 0;
489 
490  /* now sum pump rates */
491  pump_rate = 0.;
492  for( i=0; i<nN3Pump; ++i )
493  {
494  pump_rate += TauLine2[ipN3Pump[i]].Emis->pump;
495 # if 0
496  fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n",
497  i,
498  TauLine2[ipN3Pump[i]].WLAng , TauLine2[ipN3Pump[i]].pump );
499 # endif
500  }
501 
502  /* N III] N 3, N 3, 1765 multiplet */
503  /*atom_level2(&TauLines[ipT57]);*/
504  /*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
505  AtomSeqBoron(&TauLines[ipT57],
506  &TauLines[ipN3_1749],
507  &TauLines[ipN3_1747],
508  &TauLines[ipN3_1754],
509  &TauLines[ipN3_1752],
510  &TauLines[ipN3_1751],
511  0.201 , 1.088 , 0.668 , 2.044 , pump_rate,"N 3");
512  /*fprintf(ioQQQ," n4 %.3e\n", (
513  TauLines[ipN3_1749].xIntensity +
514  TauLines[ipN3_1747].xIntensity +
515  TauLines[ipN3_1754].xIntensity+
516  TauLines[ipN3_1752].xIntensity+
517  TauLines[ipN3_1751].xIntensity) / dense.xIonDense[ipNITROGEN][2] );*/
518 
519  /* N IV 1486, N 4, N 4,collisions within 3P just guess
520  * cs to ground from
521  * >>refer n4 cs Ramsbottom, C.A., Berrington, K.A., Hibbert, A., Bell, K.L. 1994,
522  * >>refercon Physica Scripta, 50, 246 */
523  if( phycon.te > 1.584e4 )
524  {
525  cs = 21.346/(phycon.te10*phycon.te10*phycon.te10*phycon.te02);
526  }
527  else
528  {
529  cs = 75.221/(phycon.sqrte/phycon.te03/phycon.te02);
530  }
531  /* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */
532  cs = MAX2(0.01,cs);
533  PutCS(cs,&TauLines[ipT1486]);
534  AtomSeqBeryllium(.9,.9,3.0,&TauLines[ipT1486],.0115);
535  embesq.em1486 = (realnum)(atoms.PopLevels[3]*0.0115*1.34e-11);
536  /*fprintf(ioQQQ," n4 %.3e\n", (TauLines[ipT1486].xIntensity + embesq.em1486 ) / dense.xIonDense[ipNITROGEN][3] );*/
537 
538  /* N IV 765, cs from
539  * >>refer n4 cs Ramsbottom, C.A., Berrington, K.A., Hibbert, A., Bell, K.L. 1994,
540  * >>refercon Physica Scripta, 50, 246 */
541  /* >>refer n4 as Flemming, J., Brage, T., Bell, K.L., Vaeck, N., Hibbert, A.,
542  * >>refercon Godefroid, M., & Froese Fischer, C., 1995, ApJ, 455, 758*/
543  cs = MIN2(4.0,1.864*phycon.te03*phycon.te03);
544  PutCS(cs,&TauLines[ipT765]);
545  atom_level2(&TauLines[ipT765]);
546 
547  /* N V 1240
548  * >>refer n5 cs Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
549  ligbar(7,&TauLines[ipT1239],&TauLines[ipT209],&cs2s2p,&cs2s3p);
550  PutCS(cs2s2p,&TauLines[ipT1239]);
551  PutCS(cs2s2p*0.5,&TauLines[ipT1243]);
552  PutCS(1.0,&TauDummy);
553  atom_level3(&TauLines[ipT1243],&TauDummy,&TauLines[ipT1239]);
554  /*fprintf(ioQQQ," n5 %.3e\n", (TauLines[ipT1243].xIntensity + TauLines[ipT1239].xIntensity ) / dense.xIonDense[ipNITROGEN][4] );*/
555 
556  /* N V 209 */
557  PutCS(cs2s3p,&TauLines[ipT209]);
558  atom_level2(&TauLines[ipT209]);
559  return;
560 }

Generated for cloudy by doxygen 1.8.1.1