cloudy trunk
Loading...
Searching...
No Matches
save_linedata.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*SaveLineData punches selected line data for all lines transferred in code */
4/*Save1LineData save data for one line */
5#include "cddefines.h"
6#include "taulines.h"
7#include "iso.h"
8#include "phycon.h"
9#include "physconst.h"
10#include "elementnames.h"
11#include "hydrogenic.h"
12#include "lines_service.h"
13#include "dense.h"
14#include "atomfeii.h"
15#include "conv.h"
16#include "lines.h"
17#include "atmdat.h"
18#include "prt.h"
19#include "h2.h"
20#include "thermal.h"
21#include "cooling.h"
22#include "save.h"
23#include "mole.h"
24
25NORETURN void SaveLineData(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
39 DEBUG_ENTRY( "SaveLineData()" );
40
41 /* routine punches out (on unit ioPUN) line data
42 * for all recombination lines, and all transitions that are transferred */
43
44 /* say what is happening so we know why we stopped */
45 fprintf( ioQQQ, " saving line data, then stopping\n" );
46
47 /* first check that all lines are turned on */
48 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
49 {
50 if( !dense.lgElmtOn[nelem] )
51 {
52 fprintf(ioQQQ," WARNING - I am saving line data but element %s is turned off.\n",
53 elementnames.chElementName[nelem]);
54 lgElemOff = true;
55 }
56 }
57 if( lgElemOff )
58 {
59 fprintf(ioQQQ,"Some elements are turned off and save line data requested.\n");
60 fprintf(ioQQQ,"Code is now designed to do save line data only with all elements on.\n");
61 fprintf(ioQQQ,"Please try again with all elements on.\n");
62 fprintf(ioQQQ,"Please try again with all elements on.\n");
64 }
65
66 /* evaluate rec coefficient at constant temperature if this is set, else
67 * use 10,000K */
68 double TeNew;
69 if( thermal.lgTemperatureConstant )
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 */
80 t_ADfA::Inst().rec_lines(phycon.te,LineSave.RecCoefCNO);
81 fprintf( ioPUN, "\n Recombination lines of C, N, O\n" );
82 fprintf( ioPUN, "#Ion\tWL(A)\tCoef\tIon\tWL(A)\tCoef\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\t%6ld\t%8.3f\t",
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.SetGasPhaseDensity( ipHYDROGEN, 1. );
101 dense.EdenHCorr = 1.;
102 EdenChange( 1. );
103
104 /* want very small neutral fractions so get mostly e- cs */
105 dense.xIonDense[ipHYDROGEN][1] = 1.e-5f;
106 findspecieslocal("H2")->den = 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 bool lgPrint = true;
136 for( i=1; i <= nLevel1; i++ )
137 {
138 /* chLineLbl generates a 1 char string from the line transfer array info -
139 * "Ne 2 128" the string is null terminated,
140 * in following printout the first 4 char is used first, followed by
141 * an integer, followed by the rest of the array*/
142 Save1LineData( TauLines[i] , ioPUN , true , lgPrint );
143 }
144
145 fprintf( ioPUN, "\n\n\n end level 1, start level 2\n" );
146 lgPrint = true;
147 for( i=0; i < nWindLine; i++ )
148 {
149 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
150 {
151 Save1LineData( TauLine2[i] , ioPUN , true , lgPrint );
152 }
153 }
154
155 fprintf( ioPUN, "\n\n\n end level 2, start inner shell UTA\n" );
156 lgPrint = true;
157 for( i=0; i < nUTA; i++ )
158 {
159 Save1LineData( UTALines[i] , ioPUN , true , lgPrint);
160 }
161
162 fprintf( ioPUN, "\n\n\n end inner shell, start H-like iso seq\n" );
163 /* h-like iso sequence */
164 /* the hydrogen like iso-sequence */
165 lgPrint = true;
166 for( nelem=0; nelem < LIMELM; nelem++ )
167 {
168 iso_collide( ipH_LIKE, nelem );
169 if( nelem < 2 || dense.lgElmtOn[nelem] )
170 {
171 /* arrays are dim'd iso_sp[ipH_LIKE][nelem].numLevels_max+1 */
172 /* keep this limit to iso.numLevels_max, instead of _local. */
173 for( ipLo=ipH1s; ipLo < iso_sp[ipH_LIKE][nelem].numLevels_max-1; ipLo++ )
174 {
175 for( ipHi=ipLo+1; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_max; ipHi++ )
176 {
177 Save1LineData( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo) , ioPUN , false,lgPrint );
178 }
179 }
180 }
181 }
182
183 fprintf( ioPUN, "\n\n\n end H-like iso seq, start He-like iso seq\n" );
184 lgPrint = true;
185 for( nelem=1; nelem < LIMELM; nelem++ )
186 {
187 if( nelem < 2 || dense.lgElmtOn[nelem] )
188 {
189 /* arrays are dim'd iso_sp[ipH_LIKE][nelem].numLevels_max+1 */
190 for( ipLo=ipHe1s1S; ipLo < iso_sp[ipHE_LIKE][nelem].numLevels_max-1; ipLo++ )
191 {
192 for( ipHi=ipLo+1; ipHi < iso_sp[ipHE_LIKE][nelem].numLevels_max; ipHi++ )
193 {
194 Save1LineData( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo) , ioPUN , false,lgPrint );
195 }
196 }
197 }
198 }
199
200 fprintf( ioPUN, "\n\n\n end he-like iso seq, start hyperfine structure lines\n" );
201 /* fine structure lines */
202 lgPrint = true;
203 for( i=0; i < nHFLines; i++ )
204 {
205 Save1LineData( HFLines[i] , ioPUN , true , lgPrint);
206 }
207
208 fprintf( ioPUN, "\n\n\n end hyperfine, start database lines\n" );
209 /* Databases: Atoms & Molecules*/
210 lgPrint = true;
211 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
212 {
213 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
214 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
215 {
216 Save1LineData( (*em).Tran() , ioPUN , true , lgPrint);
217 }
218 }
219
220 fprintf( ioPUN, "\n\n\n end database, start satellite lines\n" );
221 lgPrint = true;
222 for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
223 {
224 for( long nelem = ipISO; nelem < LIMELM; nelem++ )
225 {
226 if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
227 {
228 for( i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
229 {
230 Save1LineData( SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][i]],
231 ioPUN , true , lgPrint );
232 }
233 }
234 }
235 }
236
237 /* want very small ionized fractions so get mostly H2 cs */
238 dense.SetGasPhaseDensity( ipHYDROGEN, 1e-6f );
239 dense.EdenHCorr = 1e-6f;
240 findspecieslocal("H2")->den = 1.;
241 findspecieslocal("H2*")->den = 1.;
242 dense.xIonDense[ipHYDROGEN][1] = 1e-6;
243 EdenChange( 1e-6 );
244
245 /* H2 molecule */
246 fprintf( ioPUN, "\n\n\n" );
247 fprintf( ioPUN, " end satellite, start H2 lines\n" );
248
249 /* ioPUN unit, and option to print all possible lines - false indicates
250 * save only significant lines */
251 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
252 {
253 bool lgPopsConverged;
254 double old_val, new_val;
255 (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
256 (*diatom)->H2_Punch_line_data( ioPUN, false );
257 }
258
259 /* save FeII data if atom is currently enabled */
260 fprintf( ioPUN, "\n\n\n" );
261 fprintf( ioPUN, " end H2, start FeII lines\n" );
262
263 /* ioPUN unit, and option to print all possible lines - false indicates
264 * save only significant lines */
265 FeIIPunData( ioPUN , false );
266
267 /* ChkMonitorstring is searched for by one of the scripts in the nightly run
268 * this run will terminate with no asserts but that is the correct behavior */
269 fprintf( ioQQQ , "\n The code is left in a disturbed state after creating the SAVE LINE DATA file.\n"
270 " No calculation is actually performed, only the SAVE LINE DATA file is produced.\n"
271 " Remove the SAVE LINE DATA command to do the calculation.\n\n ChkMonitorend is ok.\n" );
272
273 /* stop when done, we have done some serious damage to the code */
275}
276
277/*Save1LineData save data for one line */
278void Save1LineData( const TransitionProxy& t , FILE * ioPUN ,
279 /* flag saying whether to give collision strength too - in multi level atoms
280 * it will be not valid without a great deal more work */
281 bool lgCS_2 ,
282 /* print header line? */
283 bool &lgPrintHead )
284{
285
286 double CritDen;
287 /* these are used for parts of the line label */
288 char chLbl[11];
289
290 DEBUG_ENTRY( "Save1LineData()" );
291
292 if( lgPrintHead )
293 fprintf( ioPUN, "#Ion\tWL\tgl\tgu\tgf\tA\tCS\tn(crt)\tdamp\n" );
294 lgPrintHead = false;
295
296 if( t.ipCont() <= 0 )
297 {
298 // not a real line, just give \n
299 return;
300 }
301
306
307 /*iWL = iWavLen( t , &chUnits , &chShift );*/
308 /* ion label, like C 1 */
309 chIonLbl(chLbl , t );
310 fprintf(ioPUN,"%s\t", chLbl );
311
312 /* this is the second piece of the line label, pick up string after start */
313
314 /* the wavelength */
315 if( strcmp( save.chConPunEnr[save.ipConPun], "labl" )== 0 )
316 {
317 prt_wl( ioPUN , t.WLAng() );
318 }
319 else
320 {
321 /* this converts energy in Rydbergs into any of the other units */
322 fprintf( ioPUN , "%.5e", AnuUnit((realnum)(t.EnergyRyd())) );
323 }
324
325 fprintf( ioPUN, "\t%3ld\t%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("\t%9.2e", t.Emis().gf()));
332
333 /* Einstein A for transition */
334 fprintf( ioPUN,PrintEfmt("\t%9.2e", t.Emis().Aul()));
335
336 /* next collision strengths, use different formats depending on size
337 * of collision strength */
338 if( t.Coll().col_str() > 100. )
339 {
340 fprintf( ioPUN, "\t%7.1f", t.Coll().col_str() );
341 }
342 else if( t.Coll().col_str() > 10. )
343 {
344 fprintf( ioPUN, "\t%7.2f", t.Coll().col_str() );
345 }
346 else if( t.Coll().col_str() > 1. )
347 {
348 fprintf( ioPUN, "\t%7.3f", t.Coll().col_str() );
349 }
350 else if( t.Coll().col_str() > .01 )
351 {
352 fprintf( ioPUN, "\t%7.4f", t.Coll().col_str() );
353 }
354 else if( t.Coll().col_str() > 0.0 )
355 {
356 fprintf( ioPUN, "\t%.3e", t.Coll().col_str() );
357 }
358 else
359 {
360 fprintf( ioPUN, "\t%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 save cs since not computed properly */
366 if( lgCS_2 && t.Coll().col_str()> 0. )
367 {
368 CritDen = t.Emis().Aul() * (*t.Hi()).g()*phycon.sqrte / (t.Coll().col_str()*COLL_CONST);
369 }
370 else
371 {
372 CritDen = 0.;
373 }
374 fprintf( ioPUN, "\t%.3e",CritDen );
375
376 // damping constant for current conditions
377 fprintf( ioPUN,PrintEfmt("\t%9.2e", t.Emis().damp() ));
378
379 fprintf( ioPUN, "\n" );
380
381 return;
382}
void FeIIPunData(FILE *ioPUN, bool lgDoAll)
FILE * ioQQQ
Definition cddefines.cpp:7
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
double AnuUnit(realnum energy)
Definition service.cpp:173
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long nWindLine
Definition cdinit.cpp:19
realnum & col_str() const
Definition collision.h:167
EmissionProxy::iterator iterator
Definition emission.h:317
realnum & damp() const
Definition emission.h:563
realnum & gf() const
Definition emission.h:513
realnum & Aul() const
Definition emission.h:613
static t_ADfA & Inst()
Definition cddefines.h:175
CollisionProxy Coll() const
Definition transition.h:424
realnum & WLAng() const
Definition transition.h:429
qList::iterator Lo() const
Definition transition.h:392
double EnergyRyd() const
Definition transition.h:83
long & ipCont() const
Definition transition.h:450
qList::iterator Hi() const
Definition transition.h:396
EmissionList::reference Emis() const
Definition transition.h:408
double den
Definition mole.h:358
void rec_lines(double t, realnum r[][471])
void EdenChange(double EdenNew)
void CoolEvaluate(double *tot)
Definition cool_eval.cpp:45
#define NORETURN
Definition cpu.h:383
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
void iso_collide(long ipISO, long nelem)
const int ipH1s
Definition iso.h:27
const int ipHE_LIKE
Definition iso.h:63
const int ipH_LIKE
Definition iso.h:62
t_LineSave LineSave
Definition lines.cpp:5
molezone * findspecieslocal(const char buf[])
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double COLL_CONST
Definition physconst.h:229
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
t_save save
Definition save.cpp:5
NORETURN void SaveLineData(FILE *ioPUN)
void Save1LineData(const TransitionProxy &t, FILE *ioPUN, bool lgCS_2, bool &lgPrintHead)
static double ** col_str
Definition species2.cpp:29
long int nSpecies
Definition taulines.cpp:21
vector< vector< TransitionList > > SatelliteLines
Definition taulines.cpp:38
TransitionList UTALines("UTALines", &AnonStates)
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
multi_arr< int, 3 > ipSatelliteLines
Definition taulines.cpp:37
long int nUTA
Definition taulines.cpp:26
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
TransitionList TauLines("TauLines", &AnonStates)
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)