cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
punch_linedata.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 /*PunchLineData punches selected line data for all lines transferred in code */
4 /*Punch1LineData punch data for one line */
5 #include "cddefines.h"
6 #include "taulines.h"
7 #include "hmi.h"
8 #include "iso.h"
9 #include "phycon.h"
10 #include "physconst.h"
11 #include "elementnames.h"
12 #include "hydrogenic.h"
13 #include "lines_service.h"
14 #include "dense.h"
15 #include "atomfeii.h"
16 #include "lines.h"
17 #include "atmdat.h"
18 #include "prt.h"
19 #include "mole_co_atom.h"
20 #include "h2.h"
21 #include "thermal.h"
22 #include "cooling.h"
23 #include "punch.h"
24 
25 NORETURN void PunchLineData(FILE * ioPUN)
26 {
27 
28  long int i,
29  j,
30  limit ,
31  nelem ,
32  ipHi ,
33  ipLo;
34 
35  const long nskip=2; /* number of emission lines per line of output */
36  double tot;
37  bool lgElemOff=false;
38  realnum a , b; /* dummy vars to pass to rotate cooling routine */
39 
40  DEBUG_ENTRY( "PunchLineData()" );
41 
42  /* routine punches out (on unit ioPUN) line data
43  * for all recombination lines, and all transitions that are transferred */
44 
45  /* say what is happening so we know why we stopped */
46  fprintf( ioQQQ, " punching line data, then stopping\n" );
47 
48  /* first check that all lines are turned on */
49  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
50  {
51  if( !dense.lgElmtOn[nelem] )
52  {
53  fprintf(ioQQQ," WARNING - I am punching line data but element %s is turned off.\n",
55  lgElemOff = true;
56  }
57  }
58  if( lgElemOff )
59  {
60  fprintf(ioQQQ,"Some elements are turned off and punch line data requested.\n");
61  fprintf(ioQQQ,"Code is now designed to do punch line data only with all elements on.\n");
62  fprintf(ioQQQ,"Please try again with all elements on.\n");
63  cdEXIT(EXIT_FAILURE);
64  }
65 
66  /* evaluate rec coefficient at constant temperature if this is set, else
67  * use 10,000K */
68  double TeNew;
70  {
71  TeNew = thermal.ConstTemp;
72  }
73  else
74  {
75  TeNew = 1e4;
76  }
77  TempChange(TeNew , false);
78 
79  /* this is set of Dima's recombination lines */
81  fprintf( ioPUN, "\n Recombination lines of C, N, O\n" );
82  fprintf( ioPUN, " Ion WL(A) Coef Ion WL(A) Coef\n" );
83  for( i=0; i<471; i+=nskip)
84  {
85  /* nskip is set to 2 above */
86  limit = MIN2(471,i+nskip);
87  fprintf( ioPUN, " " );
88  for( j=i; j < limit; j++ )
89  {
90  fprintf( ioPUN, "%2.2s%2ld%6ld%8.3f ",
91  elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][j])-1],
92  (long)(LineSave.RecCoefCNO[0][j]-LineSave.RecCoefCNO[1][j]+1.01),
93  (long)(LineSave.RecCoefCNO[2][j]+0.5),
94  log10(SDIV(LineSave.RecCoefCNO[3][j]) ) );
95  }
96  fprintf( ioPUN, " \n" );
97  }
98  fprintf( ioPUN, "\n\n" );
99 
100  dense.eden = 1.;
101  dense.gas_phase[ipHYDROGEN] = 1.;
102  dense.EdenHCorr = 1.;
103 
104  /* want very small neutral fractions so get mostly e- cs */
105  dense.xIonDense[ipHYDROGEN][1] = 1.e-5f;
106  hmi.Hmolec[ipMH2g] = 0.;
107  dense.xIonDense[ipHYDROGEN][1] = 1.;
108  for( i=1; i <= nLevel1; i++ )
109  {
110  TauLines[i].Lo->Pop = 1.;
111  }
112 
113  for( i=0; i < nWindLine; i++ )
114  {
115  TauLine2[i].Lo->Pop = 1.;
116  }
117 
118  for( i=0; i < nUTA; i++ )
119  {
120  UTALines[i].Lo->Pop = 1.;
121  }
122 
123  for( i=0; i < LIMELM; i++ )
124  {
125  for( j=0; j < LIMELM+1; j++ )
126  {
127  dense.xIonDense[i][j] = 1.;
128  }
129  }
130 
131  /* evaluate cooling, this forces evaluation of collision strengths */
132  CoolEvaluate(&tot);
133 
134  fprintf( ioPUN, " Level 1 transferred lines\n" );
135 
136  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
137 
138  for( i=1; i <= nLevel1; i++ )
139  {
140  /* chLineLbl generates a 1 char string from the line transfer array info -
141  * "Ne 2 128" the string is null terminated,
142  * in following printout the first 4 char is used first, followed by
143  * an integer, followed by the rest of the array*/
144  Punch1LineData( &TauLines[i] , ioPUN , true);
145  }
146 
147  fprintf( ioPUN, "\n\n\n" );
148  fprintf( ioPUN, " end level 1, start level 2\n" );
149  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
150  for( i=0; i < nWindLine; i++ )
151  {
152  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
153  {
154  Punch1LineData( &TauLine2[i] , ioPUN , true);
155  }
156  }
157 
158  fprintf( ioPUN, "\n\n\n" );
159  fprintf( ioPUN, " end level 2, start inner shell\n" );
160  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
161 
162  for( i=0; i < nUTA; i++ )
163  {
164  if( UTALines[i].Emis->Aul > 0. )
165  Punch1LineData( &UTALines[i] , ioPUN , true);
166  }
167 
168  fprintf( ioPUN, "\n\n\n" );
169  fprintf( ioPUN, " end inner shell, start h-like iso seq\n" );
170  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
171 
172  /* h-like iso sequence */
173  /* the hydrogen like iso-sequence */
174  for( nelem=0; nelem < LIMELM; nelem++ )
175  {
176  iso_collide( ipH_LIKE, nelem );
177  if( nelem < 2 || dense.lgElmtOn[nelem] )
178  {
179  /* arrays are dim'd iso.numLevels_max[ipH_LIKE][nelem]+1 */
180  /* keep this limit to iso.numLevels_max, instead of _local. */
181  for( ipLo=ipH1s; ipLo < iso.numLevels_max[ipH_LIKE][nelem]-1; ipLo++ )
182  {
183  for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
184  {
185  Punch1LineData( &Transitions[ipH_LIKE][nelem][ipHi][ipLo] , ioPUN , false );
186  }
187  }
188  }
189  }
190 
191  fprintf( ioPUN, "\n\n\n" );
192  fprintf( ioPUN, " end h-like iso seq, start he-like iso seq\n" );
193  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
194  for( nelem=1; nelem < LIMELM; nelem++ )
195  {
196  if( nelem < 2 || dense.lgElmtOn[nelem] )
197  {
198  /* arrays are dim'd iso.numLevels_max[ipH_LIKE][nelem]+1 */
199  for( ipLo=ipHe1s1S; ipLo < iso.numLevels_max[ipHE_LIKE][nelem]-1; ipLo++ )
200  {
201  for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
202  {
203  Punch1LineData( &Transitions[ipHE_LIKE][nelem][ipHi][ipLo] , ioPUN , false );
204  }
205  }
206  }
207  }
208 
209  fprintf( ioPUN, "\n\n\n" );
210  fprintf( ioPUN, " end he-like iso seq, start hyperfine structure lines\n" );
211  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
212  /* fine structure lines */
213  for( i=0; i < nHFLines; i++ )
214  {
215  Punch1LineData( &HFLines[i] , ioPUN , true);
216  }
217 
218  fprintf( ioPUN, "\n\n\n" );
219  fprintf( ioPUN, " end hyperfine, start database lines\n" );
220  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
221  /* Databases: Atoms & Molecules*/
222  for( i=0; i < linesAdded2; i++ )
223  {
224  Punch1LineData( atmolEmis[i].tran , ioPUN , true);
225  }
226 
227  for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
228  {
229  for( long nelem = ipISO; nelem < LIMELM; nelem++ )
230  {
231  /* must always do helium even if turned off */
232  if( nelem == ipISO || dense.lgElmtOn[nelem] )
233  {
234  for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ )
235  {
236  Punch1LineData( &SatelliteLines[ipISO][nelem][i], ioPUN , true );
237  }
238  }
239  }
240  }
241 
242  /* want very small ionized fractions so get mostly H2 cs */
243  dense.eden = 1e-6;
244  dense.gas_phase[ipHYDROGEN] = 1e-6f;
245  dense.EdenHCorr = 1e-6f;
246  dense.xIonDense[ipHYDROGEN][1] = 1.;
247  hmi.Hmolec[ipMH2g] = 1.;
248  hmi.Hmolec[ipMH2s] = 1.;
249  dense.xIonDense[ipHYDROGEN][1] = 1e-6f;
250 
251  /* H2 molecule */
252  fprintf( ioPUN, "\n\n\n" );
253  fprintf( ioPUN, " end database, start H2 lines\n" );
254  fprintf( ioPUN, "Eu Vu Ju El Vl Jl WL gl gu gf A CS n(crt)\n" );
255 
256  /* ioPUN unit, and option to print all possible lines - false indicates
257  * punch only significant lines */
258  H2_LevelPops();
259  H2_Punch_line_data( ioPUN , false );
260 
261  fprintf( ioPUN, "\n\n\n" );
262  fprintf( ioPUN, " end H2, start 12CO rotation lines\n" );
263  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
264  CO_PopsEmisCool(&C12O16Rotate, nCORotate ,1. , "12CO",&a,&b );
265  for( j=0; j< nCORotate; ++j )
266  {
267  Punch1LineData( &C12O16Rotate[j] , ioPUN , true);
268  }
269 
270  fprintf( ioPUN, "\n\n\n" );
271  fprintf( ioPUN, " end 12CO start 13CO rotation lines\n" );
272  fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" );
273  CO_PopsEmisCool(&C13O16Rotate, nCORotate , 1. ,"13CO",&a,&b );
274  for( j=0; j< nCORotate; ++j )
275  {
276  Punch1LineData( &C13O16Rotate[j] , ioPUN , true);
277  }
278 
279  /* punch FeII data if atom is currently enabled */
280  fprintf( ioPUN, "\n\n\n" );
281  fprintf( ioPUN, " end 13CO rotation lines, start FeII lines\n" );
282  fprintf( ioPUN, " Lo Hi Ion label WL gl gu gf A CS n(crt)\n" );
283 
284  /* ioPUN unit, and option to print all possible lines - false indicates
285  * punch only significant lines */
286  FeIIPunData( ioPUN , false );
287 
288  /* stop when done, we have done some serious damage to the code */
289  cdEXIT(EXIT_SUCCESS);
290 }
291 
292 /*Punch1LineData punch data for one line */
293 void Punch1LineData( transition * t , FILE * ioPUN ,
294  /* flag saying whether to give collision strength too - in multi level atoms
295  * it will be not valid without a great deal more work */
296  bool lgCS_2 )
297 {
298 
299  double CritDen;
300  /* these are used for parts of the line label */
301  char chLbl[11];
302 
303  DEBUG_ENTRY( "Punch1LineData()" );
304 
305  if( t->ipCont <= 0 )
306  {
307  return;
308  }
309 
315  /*iWL = iWavLen( t , &chUnits , &chShift );*/
316  /* ion label, like C 1 */
317  chIonLbl(chLbl , t );
318  fprintf(ioPUN,"%s\t", chLbl );
319 
320  /* this is the second piece of the line label, pick up string after start */
321 
322  /* the wavelength */
323  prt_wl(ioPUN, t->WLAng );
324 
325  fprintf( ioPUN, " %3ld%3ld",
326  /* lower and upper stat weights */
327  (long)(t->Lo->g),
328  (long)(t->Hi->g) );
329 
330  /* oscillator strength */
331  fprintf( ioPUN,PrintEfmt("%9.2e", t->Emis->gf));
332 
333  /* Einstein A for transition */
334  fprintf( ioPUN,PrintEfmt("%9.2e", t->Emis->Aul));
335 
336  /* next collision strengths, use different formats depending on size
337  * of collision strength */
338  if( t->Coll.cs > 100. )
339  {
340  fprintf( ioPUN, "%7.1f", t->Coll.cs );
341  }
342  else if( t->Coll.cs > 10. )
343  {
344  fprintf( ioPUN, "%7.2f", t->Coll.cs );
345  }
346  else if( t->Coll.cs > 1. )
347  {
348  fprintf( ioPUN, "%7.3f", t->Coll.cs );
349  }
350  else if( t->Coll.cs > .01 )
351  {
352  fprintf( ioPUN, "%7.4f", t->Coll.cs );
353  }
354  else if( t->Coll.cs > 0.0 )
355  {
356  fprintf( ioPUN, " %.3e", t->Coll.cs );
357  }
358  else
359  {
360  fprintf( ioPUN, "%7.4f", 0. );
361  }
362 
363  /* now print critical density but only if cs is positive
364  * >>chng 06 mar 24, add flag lgCS_2 - in multi-level systems do not want
365  * to punch cs since not computed properly */
366  if( lgCS_2 && t->Coll.cs> 0. )
367  {
368  CritDen = t->Emis->Aul * t->Hi->g*phycon.sqrte / (t->Coll.cs*COLL_CONST);
369  CritDen = log10(CritDen);
370  }
371  else
372  {
373  CritDen = 0.;
374  }
375  fprintf( ioPUN, "%7.3f\n",CritDen );
376  return;
377 }

Generated for cloudy by doxygen 1.8.3.1