cloudy trunk
Loading...
Searching...
No Matches
save_do.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/*SaveDo produce save output during calculation,
4 * chTime is 'MIDL' during calculation, 'LAST' at the end */
5/*SaveNewContinuum produce the 'save new continuum' output */
6/*SaveLineStuff save optical depths or source functions for all transferred lines */
7/*Save1Line called by SaveLineStuff to produce output for one line */
8/*SaveNewContinuum produce the 'save new continuum' output */
9/*SaveLineIntensity produce the 'save lines intensity' output */
10/* save h emission, for AGN3 chapter 4, routine is below */
11/*SaveResults1Line do single line of output for the save results and save line intensity commands */
12/* the number of emission lines across one line of printout */
13/*SaveSpecial generate output for the save special command */
14/*SaveResults save results from save results command */
15/*SaveResults1Line do single line of output for the save results and save line intensity commands */
16/*SaveGaunts called by save gaunts command to output gaunt factors */
17/*SaveFeII_cont return appropriate FeII banded continuum, inward, outward, or total */
18/*FindStrongestLineLabels find strongest lines contributing to point in continuum energy mesh, output in some save commands */
19#include "cddefines.h"
20#include "cddrive.h"
21#include "physconst.h"
22#include "mean.h"
23#include "taulines.h"
24#include "struc.h"
25#include "iso.h"
26#include "mole.h"
27#include "hyperfine.h"
28#include "rt.h"
29#include "lines_service.h"
30#include "doppvel.h"
31#include "dense.h"
32#include "h2.h"
33#include "magnetic.h"
34#include "hydrogenic.h"
35#include "secondaries.h"
36#include "grainvar.h"
37#include "lines.h"
38#include "dynamics.h"
39#include "colden.h"
40#include "continuum.h"
41#include "ionbal.h"
42#include "yield.h"
43#include "prt.h"
44#include "iterations.h"
45#include "heavy.h"
46#include "conv.h"
47#include "geometry.h"
48#include "called.h"
49#include "helike.h"
50#include "opacity.h"
51#include "rfield.h"
52#include "phycon.h"
53#include "timesc.h"
54#include "radius.h"
55#include "atomfeii.h"
56#include "monitor_results.h"
57#include "thermal.h"
58#include "wind.h"
59#include "hmi.h"
60#include "pressure.h"
61#include "elementnames.h"
62#include "ipoint.h"
63#include "gammas.h"
64#include "atmdat.h"
65#include "hcmap.h"
66#include "input.h"
67#include "save.h"
68#include "optimize.h"
69#include "warnings.h"
70#include "grid.h"
71#include "mole_priv.h"
72#include "atomfeii.h"
73
74// this needs C linkage since it is passed as a pointer to a C library routine
75extern "C" int wavelength_compare (const void * a, const void * b)
76{
77 LinSv a1 = *(LinSv*)a;
78 LinSv b1 = *(LinSv*)b;
79 /* comparison is b-a so we get inverse wavelength order (increasing energy order) */
80 if( b1.wavelength > a1.wavelength )
81 return 1;
82 else if( b1.wavelength < a1.wavelength )
83 return -1;
84 else
85 return 0;
86}
87
88// find strongest lines contributing to point in continuum energy mesh, output in some save commands
89STATIC void FindStrongestLineLabels( void )
90{
91 long low_index=0;
92 long high_index=0;
93 long j_min = 0;
94 double MaxFlux = 0.;
95 long ipMaxFlux = 0;
96 long j = 0;
97
98 ASSERT( LineSave.ipass==1 );
99 // copy and sort
100 memcpy(LineSvSortWL, LineSv, (unsigned)LineSave.nsum*sizeof( LinSv ));
101 qsort((void *)LineSvSortWL,(size_t)LineSave.nsum,sizeof( LinSv ),wavelength_compare);
102
103 while( rfield.anu[j_min]+0.5*rfield.widflx[j_min] < RYDLAM/LineSvSortWL[0].wavelength)
104 j_min++;
105
106 for( j=0; j<rfield.nflux; j++ )
107 {
108 if( j < j_min )
109 {
110 strcpy( rfield.chLineLabel[j], " " );
111 continue;
112 }
113
114 ASSERT( LineSvSortWL[low_index].wavelength != 0. );
115
116 while( RYDLAM/LineSvSortWL[low_index].wavelength < rfield.anu[j]-0.5*rfield.widflx[j] && low_index < LineSave.nsum-1 )
117 {
118 low_index++;
119 if( LineSvSortWL[low_index].wavelength == 0. )
120 {
121 // hit the end of real wavelengths. Pad rest of labels with spaces
122 for( long j1=j; j1<rfield.nflux; j1++ )
123 strcpy( rfield.chLineLabel[j1], " " );
124 return;
125 }
126 }
127
128 high_index = low_index;
129 ASSERT( LineSvSortWL[high_index].wavelength != 0. );
130
131 while( RYDLAM/LineSvSortWL[high_index].wavelength < rfield.anu[j]+0.5*rfield.widflx[j] && high_index < LineSave.nsum-1 )
132 {
133 high_index++;
134 if( LineSvSortWL[high_index].wavelength == 0. )
135 {
136 high_index--;
137 break;
138 }
139 }
140 // while loop found first one greater than j bin, decrement again to get back into j bin
141 high_index--;
142
143 ASSERT( LineSvSortWL[low_index].wavelength > 0. );
144 ASSERT( LineSvSortWL[high_index].wavelength > 0. );
145 ASSERT( RYDLAM/LineSvSortWL[low_index].wavelength >= rfield.anu[j]-0.5*rfield.widflx[j] );
146 ASSERT( RYDLAM/LineSvSortWL[high_index].wavelength <= rfield.anu[j]+0.5*rfield.widflx[j] );
147
148 MaxFlux = 0.;
149 ipMaxFlux = 0;
150
151 for( long k = low_index; k <= high_index; k++ )
152 {
153 if( strcmp( LineSvSortWL[k].chALab, "Coll" )==0 ||
154 strcmp( LineSvSortWL[k].chALab, "Heat" )==0 ||
155 strcmp( LineSvSortWL[k].chALab, "Pump" )==0 ||
156 strcmp( LineSvSortWL[k].chALab, "pump" )==0 ||
157 strcmp( LineSvSortWL[k].chALab, "nInu" )==0 ||
158 strcmp( LineSvSortWL[k].chALab, "nFnu" )==0 ||
159 strcmp( LineSvSortWL[k].chALab, "InwT" )==0 ||
160 strcmp( LineSvSortWL[k].chALab, "InwC" )==0 ||
161 strcmp( LineSvSortWL[k].chALab, "Inwd" )==0 ||
162 strcmp( LineSvSortWL[k].chALab, "Ca A" )==0 ||
163 strcmp( LineSvSortWL[k].chALab, "Ca B" )==0 ||
164 strcmp( LineSvSortWL[k].chALab, "Pho+" )==0 ||
165 strcmp( LineSvSortWL[k].chALab, "Pcon" )==0 ||
166 strcmp( LineSvSortWL[k].chALab, "Q(H)" )==0 ||
167 strcmp( LineSvSortWL[k].chALab, "Unit" )==0 )
168 continue;
169
170 if( LineSvSortWL[k].SumLine[0] > MaxFlux )
171 {
172 MaxFlux = LineSvSortWL[k].SumLine[0];
173 ipMaxFlux = k;
174 }
175 }
176
177 /* line label */
178 if( ipMaxFlux > 0 )
179 strcpy( rfield.chLineLabel[j], LineSvSortWL[ipMaxFlux].chALab );
180 }
181
182 return;
183}
184
185# ifdef USE_NLTE7
186//Run "Save NLTE"
187static void runNLTE(int ipPun)
188{
189 long nelem = ipNEON;
190 int nouter[11] = {0,2,3,3,3,4,3,5,5,6,0};
191 ASSERT(dense.gas_phase[nelem] != 0.);
192 // Section 1
193 fprintf( save.ipPnunit[ipPun], "data Ferland University of Kentucky Cloudy 10 11-5-2011\n");
194 fprintf( save.ipPnunit[ipPun], "case NeXY6\n");
195 fprintf( save.ipPnunit[ipPun], "code Cloudy\n");
196 fprintf( save.ipPnunit[ipPun], "atom Ne 10\n");
197 fprintf( save.ipPnunit[ipPun], "calctime %11.4e %11.4e\n",4.5,30.);
198 //Section 2
200 //Avg Charge
201 double avgchrg = 0.;
202 double m2 = 0.;
203 double m3 = 0.;
204 //double m4 = 0.;
205 for( int ion = 1; ion < nelem-1; ion++)
206 {
207 avgchrg += (dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] * ion);
208 }
209 //Central Moments of Charge Distribution
210 for( int ion = 1; ion < nelem-1; ion++)
211 {
212 m2 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW2((ion - avgchrg));
213 m3 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW3((ion - avgchrg));
214 //m4 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW4((ion - avgchrg));
215 }
216 //Specific internal energy
217 double Eint = 0.;
218 double partfun = 0.;
219 long count_nrglvls = 0;
220 for( int ion = 1; ion < nelem-1; ion++)
221 {
222 // Find ipSpecies
223 long ipSpec = 0;
224 for( int i = 0; i < nSpecies; i++)
225 {
226 ipSpec = i;
227 if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
228 break;
229 ipSpec = -1;
230 }
231
232 if( ipSpec != -1)
233 {//This requires a loop over energy levels
234 for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
235 {
236
237 Eint += dBaseStates[ipSpec][lvl].Pop()*dBaseStates[ipSpec][lvl].energy().eV();
238 partfun += dBaseStates[ipSpec][lvl].g()*exp(-1*dBaseStates[ipSpec][lvl].energy().eV()/
239 phycon.te_eV);
240 count_nrglvls += 1;
241 }
242 }
243 else
244 {
245 printf("\nNot using :\t%li\t%i\n",nelem+1,ion+1);
246 //cdEXIT(EXIT_FAILURE);
247 }
248 }
249 Eint = Eint / dense.gas_phase[nelem];
250
251
252 fprintf( save.ipPnunit[ipPun],"\nsummary_quantities\n");
253 fprintf( save.ipPnunit[ipPun],"plasma %11.4e %11.4e\n",phycon.te_eV,dense.eden);
254 fprintf( save.ipPnunit[ipPun],"time %11.4e\n",0.);
255 fprintf( save.ipPnunit[ipPun],"zbar %11.4e\n",avgchrg);
256 fprintf( save.ipPnunit[ipPun],"m2 %11.4e\n",m2);
257 fprintf( save.ipPnunit[ipPun],"m3 %11.4e\n",m3);
258 //fprintf( save.ipPnunit[ipPun],"m4 %11.3e\n",m4);
259 fprintf( save.ipPnunit[ipPun],"eint %11.4e\n",0.); //%11.3e\n",Eint);
260 fprintf( save.ipPnunit[ipPun],"deintdt %11.4e\n",0.); //%11.4e\n",thermal.dCooldT);
261 fprintf( save.ipPnunit[ipPun],"pfn %11.4e\n",partfun);
262 fprintf( save.ipPnunit[ipPun],"nmax_eff %i\n",0);
263 fprintf( save.ipPnunit[ipPun],"ploss %11.4e %11.4e %11.4e %9.3e\n",0.,0.,0.,thermal.totcol);
264 //Section 3
265 fprintf( save.ipPnunit[ipPun],"\nion_stages %li\n",nelem-2);
266
267 //Photoionization Rate
268 double PIR[nelem+1];
269 //autoionization Rate
270 double autoion[nelem+1];
271 //This requires a loop over ion stages
272 for( int ion = 1; ion < nelem-1; ion++)
273 {
274 double tempRecomb = 0.;
275 double f_aauto = 0.;
276 double f_aphoto = 0.;
277 double f_acoll = 0.;
278 double f_Scoll = 0.;
279 double f_Sphoto = 0.;
280 double f_auto = 0.;
281 autoion[ion] = 0.;
282 PIR[ion] = 0.;
283 //This deals with the fact that there are no recombinations FROM atoms
284 if( ion != 0 && ionbal.RateRecomTot[nelem][ion-1] != 0.)
285 {
286 tempRecomb = ionbal.RateRecomTot[nelem][ion-1];
287 f_aauto = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1]/tempRecomb;
288 f_aphoto = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1]/tempRecomb;
289 f_acoll = dense.eden*ionbal.CotaRate[ion-1]/tempRecomb;
290 }
291 //Photoionization rate
292 for( int i=0; i < Heavy.nsShells[nelem][ion]; i++)
293 {
294 //Add up the photoionization rates for all of the shells
295 if( ion != nelem+1)
296 {
297 PIR[ion] += ionbal.PhotoRate_Shell[nelem][ion][i][0];
298 }
299 }
300 //Compute the fractional ionizations
301 if( ion != nelem+1 && ionbal.RateIonizTot(nelem,ion) != 0.)
302 {
303 f_Scoll = ionbal.CollIonRate_Ground[nelem][ion][0]/ionbal.RateIonizTot(nelem,ion);
304 f_Sphoto = PIR[ion]/ionbal.RateIonizTot(nelem,ion);
305 f_auto = autoion[ion]/ionbal.RateIonizTot(nelem,ion);
306 }
307 //Autoionizations
308 autoion[ion] = ionbal.UTA_ionize_rate[nelem][ion] + secondaries.csupra[nelem][ion];
309 if( ion > 0)
310 {
311 autoion[ion] += 0.0; //auger[ion-1];
312 }
313 fprintf( save.ipPnunit[ipPun],"ion %li %11.4e %i\n"
314 " %11.3e %11.3e %11.3e %11.3e\n"
315 " %11.3e %11.3e %11.3e %11.3e\n\n",
316 nelem+1-ion,
317 dense.xIonDense[nelem][ion]/dense.gas_phase[nelem],
318 nouter[ion],
319 ionbal.RateIonizTot(nelem,ion),
320 f_Scoll,
321 f_Sphoto,
322 f_auto,
323 tempRecomb,
324 f_acoll,
325 f_aphoto,
326 f_aauto);
327 }
328
329 //Section 4
330 fprintf( save.ipPnunit[ipPun],"\nenergy_levels %li\n",count_nrglvls);
331 //Collisional Destruction Rate Bound-Free - Collisional Ionization
332 double GcollBF = 0.0;
333 //Photo Dest Rate Bound-Free - Photoionization
334 double GphotoBF = 0.0;
335 //Autoionization Dest Rate
336 double GautoBF = 0.0;
337 // Total Dest Rate
338 double Gtotal = 0.0;
339 //Collisional Creation Rate Bound-Free - Collisional Ionization
340 double QcollBF = 0.0;
341 //Photo Creation Rate Bound-Free - Photoionization
342 double QphotoBF = 0.0;
343 //Autoionization Creation Rate
344 double QautoBF = 0.0;
345 // Total Creation Rate
346 double Qtotal = 0.0;
347 //total pop of all levels
348 double totallevelpop = 0.;
349
350 for( int ion = 1; ion < nelem-1; ion++)
351 {
352 // Find ipSpecies
353 long ipSpec = 0;
354 for( int i = 0; i < nSpecies; i++)
355 {
356 ipSpec = i;
357 if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
358 break;
359 ipSpec = -1;
360 }
361 if( ipSpec != -1)
362 {
363 for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
364 {
365 //total level population
366 totallevelpop += dBaseStates[ipSpec][lvl].Pop();
367 }
368 }
369 }
370
371 for( int ion = 1; ion < nelem-1; ion++)
372 {
373 // Find ipSpecies
374 long ipSpec = 0;
375 for( int i = 0; i < nSpecies; i++)
376 {
377 ipSpec = i;
378 if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
379 break;
380 ipSpec = -1;
381 }
382 //Energy Level Correction Factor - Scale energy levels to ground state of atomic Neon
383 //using ionization potentials
384 double ELCF = 0.;
385 for( int i = 0; i < ion; i++ )
386 {
387 ELCF += Heavy.Valence_IP_Ryd[nelem][i]*EVRYD;
388 }
389
390 if( ipSpec != -1)
391 {
392
393 //This requires a loop over energy levels
394 for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
395 {
396
397 //Give bound-free Destruction rates for the ground state, zero otherwise
398 if( lvl == 0 && ion != nelem+1)
399 {
400 GcollBF = ionbal.CollIonRate_Ground[nelem][ion][0];
401 GphotoBF = PIR[ion];
402 GautoBF = autoion[ion];
403 }
404 else
405 {
406 GcollBF = 0.;
407 GphotoBF = 0.;
408 GautoBF = 0.;
409 }
410 //Find the Destruction total rate
411 Gtotal = dBaseStates[ipSpec][lvl].DestCollBB() +
412 dBaseStates[ipSpec][lvl].DestPhotoBB() +
413 GcollBF + GphotoBF + GautoBF;
414 //Store Dest Rate Fractions
415 double f_GcollBB,
416 f_GphotoBB,
417 f_GcollBF,
418 f_GphotoBF,
419 f_GautoBF;
420
421 if( Gtotal != 0.)
422 {
423 f_GcollBB = dBaseStates[ipSpec][lvl].DestCollBB()/Gtotal;
424 f_GphotoBB = dBaseStates[ipSpec][lvl].DestPhotoBB()/Gtotal;
425 f_GcollBF = GcollBF/Gtotal;
426 f_GphotoBF = GphotoBF/Gtotal;
427 f_GautoBF = GautoBF/Gtotal;
428 }
429 else
430 {
431 f_GcollBB = 0.;
432 f_GphotoBB = 0.;
433 f_GcollBF = 0.;
434 f_GphotoBF = 0.;
435 f_GautoBF = 0.;
436 }
437 //Give bound-free Creation rates for the ground state, zero otherwise
438 if( lvl == 0 && ion != 0 )
439 {
440 QautoBF = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1];
441 QphotoBF = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1];
442 QcollBF = dense.eden*ionbal.CotaRate[ion-1];
443 }
444 else
445 {
446 QautoBF = 0.;
447 QphotoBF = 0.;
448 QcollBF = 0.;
449 }
450
451 //Find the Creation total rate
452 Qtotal = dBaseStates[ipSpec][lvl].CreatCollBB() +
453 0.0 + QcollBF + QphotoBF + QautoBF;
454 //Get Creation Rate Fractions
455 double f_QcollBB,
456 f_QphotoBB,
457 f_QcollBF,
458 f_QphotoBF,
459 f_QautoBF;
460
461 if( Qtotal != 0.)
462 {
463 f_QcollBB = dBaseStates[ipSpec][lvl].CreatCollBB()/Qtotal;
464 f_QphotoBB = 0.0/Qtotal;
465 f_QcollBF = QcollBF/Qtotal;
466 f_QphotoBF = QphotoBF/Qtotal;
467 f_QautoBF = QautoBF/Qtotal;
468 }
469 else
470 {
471 f_QcollBB = 0.;
472 f_QphotoBB = 0.;
473 f_QcollBF = 0.;
474 f_QphotoBF = 0.;
475 f_QautoBF = 0.;
476 }
477
478 //Section 4 print
479 fprintf( save.ipPnunit[ipPun],"elev %li %3i %11.4e %11.6e %11.4e\n"
480 " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n"
481 " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n"
482 " %i %i %i\n\n",
483 nelem+1-ion,
484 lvl+1,
485 dBaseStates[ipSpec][lvl].g(),
486 dBaseStates[ipSpec][lvl].energy().eV()+ELCF,
487 dBaseStates[ipSpec][lvl].Pop()/totallevelpop,
488 Gtotal*dBaseStates[ipSpec][lvl].Pop(),
489 f_GcollBB,
490 f_GphotoBB,
491 f_GcollBF,
492 f_GphotoBF,
493 f_GautoBF,
494 Qtotal,
495 f_QcollBB,
496 f_QphotoBB,
497 f_QcollBF,
498 f_QphotoBF,
499 f_QautoBF,
500 0,
501 0,
502 nouter[ion]);
503 }
504 }
505 else
506 printf("\nNot using :\t%li\t%i for levels\n",nelem+1,ion+1);
507 }
508 return;
509}
510# endif
511
512
513// implements the absorption option on the
514// set save line width command
515inline realnum PrettyTransmission(long j, realnum transmission)
516{
517 if( save.ResolutionAbs < realnum(0.) )
518 // option to conserve energy
519 return transmission;
520 else
521 {
522 realnum corr = save.ResolutionAbs*rfield.widflx[j]/rfield.anu[j];
523 return realnum(max(0., 1. - (1.-transmission)*corr ));
524 }
525}
526
527/*SaveResults1Line do single line of output for the save results and save line intensity commands */
528/* the number of emission lines across one line of printout */
529STATIC void SaveResults1Line(
530 FILE * ioPUN,
531 /* 4 char + null string */
532 const char *chLab,
533 realnum wl,
534 double xInten,
535 const char *chFunction);/* null terminated string saying what to do */
536
537/*SaveGaunts called by save gaunts command to output gaunt factors */
538STATIC void SaveGaunts(FILE* ioPUN);
539
540/*SaveResults save results from save results command */
541/*SaveResults1Line do single line of output for the save results and save line intensity commands */
542STATIC void SaveResults(FILE* ioPUN);
543
544STATIC void SaveLineStuff(
545 FILE * ioPUN,
546 const char *chJob ,
547 realnum xLimit);
548
549/* save h emission, for chapter 4, routine is below */
550STATIC void AGN_Hemis(FILE *ioPUN );
551
552/*SaveNewContinuum produce the 'save new continuum' output */
553STATIC void SaveNewContinuum(FILE * ioPUN );
554
555/*SaveLineIntensity produce the 'save lines intensity' output */
556STATIC void SaveLineIntensity(FILE * ioPUN , long int ipPun, realnum Threshold);
557
558char *chDummy;
559
560/*SaveFeII_cont return appropriate FeII banded continuum, inward, outward, or total */
561STATIC realnum SaveFeII_cont( long int ipCont , long ipFeII_Cont_type )
562{
563
564 DEBUG_ENTRY( "SaveFeII_cont()" );
565
566 // [0]=wavelength, [1] inward; [2] outward, 3=> total
567 if( ipFeII_Cont_type == 3 )
568 return( FeII_Cont[ipCont][1]+FeII_Cont[ipCont][2] );
569 else
570 return(FeII_Cont[ipCont][ipFeII_Cont_type]);
571}
572
573void SaveDo(
574 /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
575 const char *chTime)
576{
577 bool lgTlkSav;
578 long int
579 ipPun,
580 i,
581 i1,
582 ion,
583 ipConMax,
584 ipHi,
585 ipLinMax,
586 ipLo,
587 ips,
588 j,
589 jj,
590 limit,
591 nelem;
592 double ConMax,
593 RateInter,
594 av,
595 conem,
596 eps,
597 flxatt,
598 flxcor,
599 flxin,
600 flxref,
601 flxtrn,
602 fout,
603 fref,
604 fsum,
605 opConSum,
606 opLinSum,
607 stage,
608 sum,
609 texc,
610 xLinMax;
611
612 DEBUG_ENTRY( "SaveDo()" );
613
614 /*
615 * the "last" option on save command, to save on last iteration,
616 * is parsed at the top of the loop in only one place.
617 * no further action is needed at all for save last to work
618 * ok throughout this routine
619 */
620
621 /*
622 * each branch can have a test whether chTime is or is not "LAST"
623 *
624 * if( strcmp(chTime,"LAST") == 0 ) <== print after iteration is complete
625 *
626 * if "LAST" then this is last call to routine after iteration complete
627 * save only if "LAST" when results at end of iteration are needed
628 *
629 * if( strcmp(chTime,"LAST") != 0 ) <== print for every zone
630 *
631 * test for .not."LAST" is for every zone result, where you do not
632 * want to save last zone twice
633 */
634
635 /* return if no save to do */
636 if( save.nsave < 1 )
637 {
638 return;
639 }
640
641 /* during a grid calculation this routine saves grid points after
642 * cloudy is called. we may output it below */
643 if( grid.lgGrid )
644 {
645 if( chTime[0]=='L' )
646 GridGatherInCloudy();
647 }
648
649 // sort line labels if this is last call, this avoids multiple calls if several
650 // output options need sorted labels and is safer since labels will be sorted in
651 // case new code is added that reports the strong lines. The disadvantage is that
652 // we sort even if the labels are not used
653 if( strcmp(chTime,"LAST") == 0 )
654 {
655 // sort emission line intensities so strongest lines are reported
656 FindStrongestLineLabels();
657 }
658
659 for( ipPun=0; ipPun < save.nsave; ipPun++ )
660 {
661 if( save.lgPunHeader[ipPun] && !nMatch(save.chHeader[ipPun],save.chNONSENSE) )
662 {
663 fprintf( save.ipPnunit[ipPun], "%s", save.chHeader[ipPun] );
664 save.lgPunHeader[ipPun] = false;
665 }
666
667 /* this global variable to remember where in the save stack we are */
668 save.ipConPun = ipPun;
669
670 /* used to identify case where no key found */
671 bool lgNoHitFirstBranch = false;
672
673 /* iterations.lgLastIt is true if this is last iteration
674 * lgPunLstIter set true if 'last' key occurred on save command
675 * normally is false. This will skip saving if last set and
676 * this is not last iteration */
677 /* IMPORTANT: there is a second, identical if-statement halfway
678 * down this routine. Any changes here should be copied there! */
679 if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
680 // if the sim aborted, make sure that the punch output is still done.
681 // for MIDL output this is not very useful as all previous zones from
682 // the iteration where the abort occured will be missing, but for LAST
683 // output it is important to print at least placeholders in grid runs.
684 ( lgAbort && strcmp(chTime,"LAST") == 0 ) ||
685 ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) )
686 {
687
688 if( strcmp(save.chSave[ipPun],"ABUN") == 0 )
689 {
690 /* save abundances vs depth */
691 if( strcmp(chTime,"LAST") != 0 )
692 {
693 fprintf( save.ipPnunit[ipPun], "%.2f",
694 log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) );
695 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
696 {
697 /* >>chng 05 feb 03, protect against non-positive abundances,
698 * bug caught by Marcelo Castellanos */
699 fprintf( save.ipPnunit[ipPun], "\t%.2f",
700 log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) );
701 }
702 fprintf( save.ipPnunit[ipPun], "\n" );
703 }
704 }
705
706 else if( strcmp(save.chSave[ipPun],"21CM") == 0 )
707 {
708 /* save information about 21 cm line */
709 if( strcmp(chTime,"LAST") != 0 )
710 {
711 fprintf( save.ipPnunit[ipPun],
712 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
713 /* depth, cm */
714 radius.depth_mid_zone,
715 hyperfine.Tspin21cm ,
716 phycon.te ,
717 /* temperature from Lya - 21 cm optical depth ratio */
718 3.84e-7* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon() /
719 SDIV( HFLines[0].Emis().TauCon() ),
720 /*TexcLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),*/
721 (*HFLines[0].Lo()).Pop() ,
722 (*HFLines[0].Hi()).Pop() ,
723 OccupationNumberLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),
724 HFLines[0].Emis().TauCon() ,
725 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon(),
726 HFLines[0].Emis().PopOpc(),
727 /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */
728 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop())/SDIV( hyperfine.Tspin21cm),
729 /* why was above multiplied by this following term? */
730 /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/
731 HFLines[0].Emis().TauIn(),
732 TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ) ,
733 colden.H0_ov_Tspin,
734 /*>>chng 27 mar, GS, integrated 21cm spin temperature*/
735 colden.H0_21cm_lower,
736 colden.H0_21cm_upper,
737 -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower)
738 );
739 }
740 }
741
742 else if( strcmp(save.chSave[ipPun],"AGES") == 0 )
743 {
744 /* save timescales vs depth */
745 if( strcmp(chTime,"LAST") != 0 )
746 {
747 int ipCO, ipOH;
748 ipCO = findspecies("CO")->index;
749 ipOH = findspecies("OH")->index;
750 fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
751 /* depth, cm */
752 radius.depth_mid_zone,
753 /* cooling timescale */
754 dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot,
755 /* H2 destruction timescale */
756 timesc.time_H2_Dest_here,
757 /* CO destruction timescale */
758 1./SDIV((ipCO != -1) ? mole.species[ipCO].snk : 0.),
759 /* OH destruction timescale */
760 1./SDIV((ipOH != -1) ? mole.species[ipOH].snk : 0.),
761 /* H recombination timescale */
762 1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) );
763 }
764 }
765
766 else if( strcmp(save.chSave[ipPun]," AGN") == 0 )
767 {
768 if( strcmp(chTime,"LAST") == 0 )
769 {
770 if( strcmp( save.chSaveArgs[ipPun], "HECS" ) == 0 )
771 {
772 /* this routine is in helike.c */
773 AGN_He1_CS(save.ipPnunit[ipPun]);
774 }
775 if( strcmp( save.chSaveArgs[ipPun], "HEMI" ) == 0 )
776 {
777 /* save h emiss, for chapter 4, routine is below */
778 AGN_Hemis(save.ipPnunit[ipPun]);
779 }
780 else
781 {
782 fprintf( ioQQQ, " SaveDo does not recognize flag %4.4s set for AGN save. This is impossible.\n",
783 save.chSave[ipPun] );
784 ShowMe();
785 cdEXIT(EXIT_FAILURE);
786 }
787 }
788 }
789
790 else if( strcmp(save.chSave[ipPun],"MONI") == 0 )
791 {
792 if( strcmp(chTime,"LAST") == 0 )
793 {
794 /* save the monitor output */
795 lgCheckMonitors(save.ipPnunit[ipPun]);
796 }
797 }
798
799 else if( strcmp(save.chSave[ipPun],"AVER") == 0 )
800 {
801 if( strcmp(chTime,"LAST") == 0 )
802 {
803 /* save the averages output */
804 save_average( ipPun );
805 }
806 }
807
808 else if( strncmp(save.chSave[ipPun],"CHA",3) == 0 )
809 {
810 if( strcmp(chTime,"LAST") == 0 )
811 {
812 /* one of the charge transfer options, all in chargtran.c */
813 ChargTranPun( save.ipPnunit[ipPun] , save.chSave[ipPun] );
814 }
815 }
816
817 else if( strcmp( save.chSave[ipPun],"CHIA") == 0)
818 {
819 static bool lgRunOnce = true;
820 if( lgRunOnce )
821 {
822 lgRunOnce = false;
823 // save chianti collision data in physical units
824 int ipLo = 0;
825 int ipHi = 0;
826 double fupsilon = 0.;
827 double initTemp = 4.0;
828 double finalTemp = 9.1;
829 double stepTemp = 0.2;
830 fprintf( save.ipPnunit[ipPun],"Species\tLo\tHi\tWLAng");
831 for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
832 {
833 fprintf( save.ipPnunit[ipPun],"\t%2.1f",logtemp);
834 }
835 fprintf( save.ipPnunit[ipPun],"\n");
836 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
837 {
838 for( EmissionList::iterator tr=dBaseTrans[ipSpecies].Emis().begin();
839 tr != dBaseTrans[ipSpecies].Emis().end(); ++tr)
840 {
841 if( dBaseTrans[ipSpecies].chLabel() == "Chianti" )
842 {
843 ipLo = tr->Tran().ipLo();
844 ipHi = tr->Tran().ipHi();
845 fprintf( save.ipPnunit[ipPun],"%s\t%i\t%i\t",
846 dBaseSpecies[ipSpecies].chLabel,ipLo+1,ipHi+1);
847 fprintf( save.ipPnunit[ipPun],"%5.3e",tr->Tran().WLAng());
848 for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
849 {
850 fupsilon = CHIANTI_Upsilon(ipSpecies, ipELECTRON, ipHi, ipLo, pow(10.,logtemp));
851 fprintf( save.ipPnunit[ipPun],"\t%5.3e",fupsilon);
852 }
853 fprintf( save.ipPnunit[ipPun],"\n");
854 }
855 }
856 }
857 }
858 }
859
860 else if( strcmp(save.chSave[ipPun],"COOL") == 0 ||
861 strcmp(save.chSave[ipPun],"EACH") == 0 )
862 {
863 /* save cooling, routine in file of same name */
864 if( strcmp(chTime,"LAST") != 0 )
865 CoolSave(save.ipPnunit[ipPun], save.chSave[ipPun]);
866 }
867
868 else if( strcmp(save.chSave[ipPun],"DOMI") == 0 )
869 {
870 /* save dominant rates */
871 if( strcmp(chTime,"LAST") != 0 )
872 {
873 molecule *debug_species = findspecies( save.chSpeciesDominantRates[ipPun] );
874 if (debug_species == null_mole)
875 {
876 fprintf( ioQQQ,"Error in SAVE DOMINANT RATES, species %s not found\n",
877 save.chSpeciesDominantRates[ipPun]);
878 }
879 else
880 {
881 fprintf( save.ipPnunit[ipPun],
882 "%e\t%e\t", radius.depth_mid_zone, mole.species[ debug_species->index ].column );
883 mole_dominant_rates( debug_species, save.ipPnunit[ipPun] );
884 }
885 }
886 }
887
888 else if( strncmp(save.chSave[ipPun],"DYN" , 3) == 0 )
889 {
890 /* save dynamics xxx, information about dynamical solutions */
891 if( strcmp(chTime,"LAST") != 0 )
892 DynaSave(save.ipPnunit[ipPun] ,save.chSave[ipPun][3] );
893 }
894
895 else if( strcmp(save.chSave[ipPun],"ENTH") == 0 )
896 {
897 if( strcmp(chTime,"LAST") != 0 )
898 fprintf( save.ipPnunit[ipPun],
899 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
900 radius.depth_mid_zone,
901 phycon.EnthalpyDensity,
902 phycon.EnergyExcitation,
903 phycon.EnergyIonization,
904 phycon.EnergyBinding ,
905 0.5*POW2(wind.windv)*dense.xMassDensity , /* KE */
906 5./2.*pressure.PresGasCurr , /* thermal plus PdV work */
907 magnetic.EnthalpyDensity); /* magnetic terms */
908 }
909
910 else if( strcmp(save.chSave[ipPun],"COLU") == 0 )
911 {
912 /* save column densities */
913 if( strcmp(chTime,"LAST") == 0 )
914 {
915 PrtColumns(save.ipPnunit[ipPun],"TABLE",ipPun);
916 }
917 }
918
919 else if( strcmp(save.chSave[ipPun],"COLS") == 0 )
920 {
921 /* save some column densities */
922 if( strcmp(chTime,"LAST") == 0 )
923 {
924 /*TODO unify following with PrtColumns above */
925 save_colden( save.ipPnunit[ipPun] );
926 }
927 }
928
929 else if( strcmp(save.chSave[ipPun],"COMP") == 0 )
930 {
931 /* Compton energy exchange coefficients */
932 if( strcmp(chTime,"LAST") != 0 )
933 {
934 for( jj=0; jj<rfield.nflux; jj = jj + save.ncSaveSkip)
935 {
936 fprintf( save.ipPnunit[ipPun], "%10.2e%10.2e%10.2e\n",
937 rfield.anu[jj], rfield.comup[jj]/rfield.widflx[jj],
938 rfield.comdn[jj]/rfield.widflx[jj] );
939 }
940 }
941 }
942
943 /* save continuum commands */
944 else if( strcmp(save.chSave[ipPun],"CON ") == 0 )
945 {
946 /* this is the usual "save continuum" case */
947 /* >>chng 06 apr 03, add every option to do every zone */
948 /* if save every is set then nSaveEveryZone must be positive
949 * was init to -1 */
950 bool lgPrintThis =false;
951 if( save.lgSaveEveryZone[ipPun] )
952 {
953 /* this branch, every option is on line so want to print every n zone */
954 if( strcmp(chTime,"LAST") != 0 )
955 {
956 /* not last zone - print first and intermediate cases */
957 if( nzone==1 )
958 {
959 lgPrintThis = true;
960 }
961 else if( nzone%save.nSaveEveryZone[ipPun]==0 )
962 {
963 lgPrintThis = true;
964 }
965 }
966 else
967 {
968 /* this is last zone, print only if did not trip on above */
969 if( nzone!=1 && nzone%save.nSaveEveryZone[ipPun]!=0 )
970 {
971 lgPrintThis = true;
972 }
973 }
974 }
975 else
976 {
977 /* this branch, not "every", so only print the last zone */
978 if( strcmp(chTime,"LAST") == 0 )
979 lgPrintThis = true;
980 }
981 ASSERT( !save.lgSaveEveryZone[ipPun] || save.nSaveEveryZone[ipPun]>0 );
982 if( lgPrintThis )
983 {
984 if( save.lgSaveEveryZone[ipPun] && nzone!=1)
985 fprintf( save.ipPnunit[ipPun], "%s\n",
986 save.chHashString );
987
988 /* option to also print same arrays but for cumulative arrays */
989 int nEmType = (int)save.punarg[ipPun][0];
990 ASSERT( nEmType==0 || nEmType==1 );
991
992 const realnum *trans_coef_total=rfield.getCoarseTransCoef();
993 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
994 {
995 /* four continua predicted here;
996 * incident, attenuated incident, emitted,
997 * then attenuated incident + emitted, last reflected
998 * reflected continuum is stored relative to inner radius
999 * others are stored for this radius */
1000
1001 /* NB this code also used in save emitted,
1002 * transmitted continuum commands */
1003
1004 /* the incident continuum */
1005 flxin = ((double)rfield.flux_total_incident[nEmType][j])
1006 *rfield.anu2[j]*
1007 EN1RYD/rfield.widflx[j];
1008
1009 // a value < 0. indicates that energy should be conserved
1010 realnum resolution = ( save.Resolution < realnum(0.) ) ?
1011 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1012
1013 /* the reflected continuum */
1014 flxref = (rfield.anu2[j]*((double)rfield.ConRefIncid[nEmType][j]+rfield.ConEmitReflec[nEmType][j])/rfield.widflx[j] +
1015 rfield.anu[j]*resolution*rfield.reflin[nEmType][j])*EN1RYD;
1016
1017 /* the attenuated incident continuum */
1018 flxatt = ((double)rfield.flux[nEmType][j])*rfield.anu2[j]*EN1RYD/
1019 rfield.widflx[j]*radius.r1r0sq *
1020 PrettyTransmission( j, trans_coef_total[j] );
1021
1022 /* the outward emitted continuum */
1023 conem = ((double)rfield.ConEmitOut[nEmType][j]/
1024 rfield.widflx[j]*rfield.anu2[j] + resolution*
1025 rfield.outlin[nEmType][j]*rfield.anu[j])*radius.r1r0sq*
1026 EN1RYD*geometry.covgeo;
1027
1028 /* sum of emitted and transmitted continua */
1029 flxtrn = conem + flxatt;
1030
1031 /* photon energy in appropriate energy or wavelength units */
1032 fprintf( save.ipPnunit[ipPun],"%.5e\t", AnuUnit(rfield.AnuOrg[j]) );
1033 /* incident continuum */
1034 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxin );
1035 /* trans cont */
1036 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxatt );
1037 /* DiffOut cont */
1038 fprintf( save.ipPnunit[ipPun],"%.3e\t", conem );
1039 /* net trans cont */
1040 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxtrn );
1041 /* reflected cont */
1042 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxref );
1043 /* total cont */
1044 fprintf( save.ipPnunit[ipPun],"%.3e\t", flxref + flxtrn );
1045
1046 /* reflected lines */
1047 fprintf( save.ipPnunit[ipPun],"%.3e\t", rfield.anu[j]*resolution*rfield.reflin[nEmType][j]*EN1RYD );
1048
1049 /* outward lines */
1050 fprintf( save.ipPnunit[ipPun],"%.3e\t", rfield.anu[j]*resolution*rfield.outlin[nEmType][j]*EN1RYD*
1051 radius.r1r0sq*geometry.covgeo );
1052
1053 fprintf( save.ipPnunit[ipPun], "%4.4s\t%4.4s\t",
1054 /* line label */
1055 rfield.chLineLabel[j] ,
1056 /* cont label*/
1057 rfield.chContLabel[j] );
1058 /* number of lines within that cell over cell width
1059 * save raw continuum has number by itself */
1060 fprintf( save.ipPnunit[ipPun], "%.2f\n", rfield.line_count[j]/rfield.widflx[j]*rfield.anu[j] );
1061 }
1062 }
1063 }
1064
1065 /* this is the save spectrum command,
1066 * the new continuum command that will replace the previous one */
1067 else if( strcmp(save.chSave[ipPun],"CONN") == 0 )
1068 {
1069 if( strcmp(chTime,"LAST") == 0 )
1070 SaveNewContinuum( save.ipPnunit[ipPun]);
1071 }
1072
1073 else if( strcmp(save.chSave[ipPun],"CONC") == 0 )
1074 {
1075 /* save incident continuum */
1076 /* set pointer for possible change in units of energy in continuum
1077 * AnuUnit will give anu in whatever units were set with save units */
1078 if( strcmp(chTime,"LAST") == 0 )
1079 {
1080 /* incident continuum */
1081 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1082 {
1083 flxin = rfield.flux_total_incident[0][j]*rfield.anu2[j]*
1084 EN1RYD/rfield.widflx[j];
1085 /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
1086 fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\n",
1087 AnuUnit(rfield.AnuOrg[j]), flxin );
1088 }
1089 }
1090 }
1091
1092 else if( strcmp(save.chSave[ipPun],"CONG") == 0 )
1093 {
1094 /* save emitted grain continuum in optically thin limit */
1095 if( strcmp(chTime,"LAST") == 0 )
1096 {
1097 /* GrainMakeDiffuse broke out emission into types
1098 * according to matType */
1099 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1100 {
1101 fsum = gv.GraphiteEmission[j]*rfield.anu2[j]*
1102 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
1103
1104 fout = gv.SilicateEmission[j]*rfield.anu2[j]*
1105 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
1106
1107 /* anu is .5e format to resolve energy mesh near 1 Ryd
1108 * AnuUnit gives anu in whatever units were set with units option */
1109 fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n",
1110 AnuUnit(rfield.AnuOrg[j]) , fsum , fout ,
1111 /* total emission per unit volume defined in GrainMakeDiffuse
1112 * used in RT_diffuse to form total grain emission */
1113 gv.GrainEmission[j]*rfield.anu2[j]*
1114 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo );
1115 }
1116 }
1117 }
1118
1119 else if( strcmp(save.chSave[ipPun],"CONR") == 0 )
1120 {
1121 /* save reflected continuum */
1122 /* set pointer for possible change in units of energy in continuum
1123 * AnuUnit will give anu in whatever units were set with save units */
1124 if( strcmp(chTime,"LAST") == 0 )
1125 {
1126 if( geometry.lgSphere )
1127 {
1128 fprintf( save.ipPnunit[ipPun], " Reflected continuum not predicted when SPHERE is set.\n" );
1129 fprintf( ioQQQ ,
1130 "\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" );
1131 cdEXIT(EXIT_FAILURE);
1132 }
1133
1134 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1135 {
1136 // a value < 0. indicates that energy should be conserved
1137 realnum resolution = ( save.Resolution < realnum(0.) ) ?
1138 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1139
1140 /* reflected continuum */
1141 flxref = rfield.anu2[j]*((double)rfield.ConRefIncid[0][j]+rfield.ConEmitReflec[0][j])/
1142 rfield.widflx[j]*EN1RYD;
1143 /* reflected lines */
1144 fref = rfield.anu[j]*resolution*rfield.reflin[0][j]*EN1RYD;
1145 /* ratio of reflected to incident continuum, the albedo */
1146 if( rfield.flux_total_incident[0][j] > 1e-25 )
1147 {
1148 av = rfield.ConRefIncid[0][j]/rfield.flux_total_incident[0][j];
1149 }
1150 else
1151 {
1152 av = 0.;
1153 }
1154 /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
1155 fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4s\n",
1156 AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref,
1157 av, rfield.chContLabel[j] );
1158 }
1159 }
1160 }
1161
1162 else if( strcmp(save.chSave[ipPun],"CNVE") == 0 )
1163 {
1164 /* the save convergence error command */
1165 if( strcmp(chTime,"LAST") != 0 )
1166 {
1167 fprintf( save.ipPnunit[ipPun],
1168 "%.5e\t%li\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n",
1169 radius.depth_mid_zone,
1170 conv.nPres2Ioniz,
1171 pressure.PresTotlCurr,
1172 pressure.PresTotlError*100.,
1173 dense.EdenTrue,
1174 dense.eden,
1175 (dense.EdenTrue - dense.eden)*100./dense.EdenTrue,
1176 thermal.htot,
1177 thermal.ctot,
1178 (thermal.htot - thermal.ctot)*100./thermal.htot );
1179 }
1180 }
1181
1182 else if( strcmp(save.chSave[ipPun],"CONB") == 0 )
1183 {
1184 /* save continuum bins binning */
1185 /* set pointer for possible change in units of energy in continuum
1186 * AnuUnit will give anu in whatever units were set with save units */
1187 if( strcmp(chTime,"LAST") != 0 )
1188 {
1189 for( j=0; j<rfield.nupper; j = j + save.ncSaveSkip)
1190 {
1191 fprintf( save.ipPnunit[ipPun], "%14.5e %14.5e %14.5e\n",
1192 AnuUnit(rfield.AnuOrg[j]), rfield.anu[j], rfield.widflx[j] );
1193 }
1194 }
1195 }
1196
1197 else if( strcmp(save.chSave[ipPun],"COND") == 0 )
1198 {
1199 ASSERT( fp_equal( save.punarg[ipPun][0] , (realnum)2. ) ||
1200 fp_equal( save.punarg[ipPun][0] , (realnum)1.) );
1201
1202 /* save diffuse continuum the local line and continuous emission */
1203 if( strcmp(chTime,"LAST") == 0 &&
1204 fp_equal(save.punarg[ipPun][0] , (realnum)1.) )
1205 {
1206 /* this option to save diffuse emission for last zone
1207 * as series of columns */
1208 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1209 {
1210 // a value < 0. indicates that energy should be conserved
1211 realnum resolution = ( save.Resolution < realnum(0.) ) ?
1212 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1213
1214 double EmisLin = resolution*EN1RYD*
1215 rfield.DiffuseLineEmission[j]*rfield.anu[j];
1216 double EmisCon = rfield.ConEmitLocal[nzone][j]*
1217 rfield.anu2[j]*EN1RYD/rfield.widflx[j];
1218 fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\t%.5e\n",
1219 AnuUnit(rfield.AnuOrg[j]),
1220 EmisCon ,
1221 EmisLin ,
1222 EmisCon+EmisLin);
1223 }
1224 }
1225 else if( strcmp(chTime,"LAST") != 0 &&
1226 fp_equal(save.punarg[ipPun][0] , (realnum)2.) )
1227 {
1228 /* diffuse emission per zone, with each a long row */
1229 static bool lgMustPrintHeader=true;
1230 if( lgMustPrintHeader )
1231 {
1232 lgMustPrintHeader = false;
1233 fprintf( save.ipPnunit[ipPun], "%.5e",
1234 AnuUnit(rfield.AnuOrg[0]) );
1235 for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
1236 {
1237 fprintf( save.ipPnunit[ipPun], "\t%.5e",
1238 AnuUnit(rfield.AnuOrg[j]) );
1239 }
1240 fprintf( save.ipPnunit[ipPun], "\n" );
1241 }
1242 // a value < 0. indicates that energy should be conserved
1243 realnum resolution = ( save.Resolution < realnum(0.) ) ?
1244 rfield.anu[0]/rfield.widflx[0] : save.Resolution;
1245 double EmisLin = resolution*EN1RYD*
1246 rfield.DiffuseLineEmission[0]*rfield.anu[0];
1247 double EmisCon = rfield.ConEmitLocal[nzone][0]*
1248 rfield.anu2[0]*EN1RYD/rfield.widflx[0];
1249 fprintf( save.ipPnunit[ipPun], "%.5e",
1250 EmisCon+EmisLin);
1251 for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
1252 {
1253 // a value < 0. indicates that energy should be conserved
1254 resolution = ( save.Resolution < realnum(0.) ) ?
1255 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1256 double EmisLin = resolution*EN1RYD*
1257 rfield.DiffuseLineEmission[j]*rfield.anu[j];
1258 double EmisCon = rfield.ConEmitLocal[nzone][j]*
1259 rfield.anu2[j]*EN1RYD/rfield.widflx[j];
1260 fprintf( save.ipPnunit[ipPun], "\t%.5e",
1261 EmisCon+EmisLin);
1262 }
1263 fprintf( save.ipPnunit[ipPun], "\n" );
1264 }
1265 }
1266
1267 else if( strcmp(save.chSave[ipPun],"CONE") == 0 )
1268 {
1269 /* save emitted continuum */
1270 /* set pointer for possible change in units of energy in continuum
1271 * AnuUnit will give anu in whatever units were set with save units */
1272 if( strcmp(chTime,"LAST") == 0 )
1273 {
1274 /* save emitted continuum */
1275 for( j=0; j<rfield.nflux;j = j +save.ncSaveSkip)
1276 {
1277 // a value < 0. indicates that energy should be conserved
1278 realnum resolution = ( save.Resolution < realnum(0.) ) ?
1279 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1280
1281 /* this is the reflected component */
1282 flxref = (rfield.anu2[j]*(rfield.ConRefIncid[0][j]+rfield.ConEmitReflec[0][j])/
1283 rfield.widflx[j] + rfield.anu[j]*resolution*
1284 rfield.reflin[0][j])*EN1RYD;
1285
1286 /* this is the total emission in the outward direction */
1287 conem = (rfield.ConEmitOut[0][j])/rfield.widflx[j]*rfield.anu2[j] +
1288 resolution*rfield.outlin[0][j]*rfield.anu[j];
1289
1290 conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;
1291
1292 /* output: photon energy, reflected, outward, total emission
1293 * >>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */
1294 fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n",
1295 AnuUnit(rfield.AnuOrg[j]),
1296 flxref,
1297 conem,
1298 flxref + conem,
1299 rfield.chLineLabel[j],
1300 rfield.chContLabel[j]
1301 );
1302 }
1303 }
1304 }
1305
1306 /* save fine continuum command */
1307 else if( strcmp(save.chSave[ipPun],"CONf") == 0 )
1308 {
1309 if( strcmp(chTime,"LAST") == 0 )
1310 {
1311 long nu_hi , nskip;
1312 if( save.punarg[ipPun][0] > 0. )
1313 /* possible lower bounds to energy range -
1314 * 0 if not set with range option*/
1315 j = ipFineCont( save.punarg[ipPun][0] );
1316 else
1317 j = 0;
1318
1319 /* upper limit set with range option */
1320 if( save.punarg[ipPun][1]> 0. )
1321 nu_hi = ipFineCont( save.punarg[ipPun][1]);
1322 else
1323 nu_hi = rfield.nfine;
1324
1325 /* number of cells to bring together, default is 10 */
1326 nskip = (long)save.punarg[ipPun][2];
1327
1328 do
1329 {
1330 realnum sum1 = rfield.fine_opt_depth[j];
1331 realnum xnu = rfield.fine_anu[j];
1332 for( jj=1; jj<nskip; ++jj )
1333 {
1334 xnu += rfield.fine_anu[j+jj];
1335 sum1 += rfield.fine_opt_depth[j+jj];
1336 }
1337 fprintf( save.ipPnunit[ipPun],
1338 "%.6e\t%.3e\n",
1339 AnuUnit(xnu/nskip),
1340 sexp(sum1/nskip) );
1341 j += nskip;
1342 } while( j < nu_hi );
1343 }
1344 }
1345
1346 else if( strcmp(save.chSave[ipPun],"CONi") == 0 )
1347 {
1348 /* save continuum interactions */
1349 /* set pointer for possible change in units of energy in continuum
1350 * AnuUnit will give anu in whatever units were set with save units */
1351
1352 /* continuum interactions */
1353 if( strcmp(chTime,"LAST") != 0 )
1354 {
1355 /* this is option to set lowest energy */
1356 if( save.punarg[ipPun][0] <= 0. )
1357 {
1358 i1 = 1;
1359 }
1360 else if( save.punarg[ipPun][0] < 100. )
1361 {
1362 i1 = ipoint(save.punarg[ipPun][0]);
1363 }
1364 else
1365 {
1366 i1 = (long int)save.punarg[ipPun][0];
1367 }
1368
1369 fref = 0.;
1370 fout = 0.;
1371 fsum = 0.;
1372 sum = 0.;
1373 flxin = 0.;
1374
1375 for( j=i1-1; j < rfield.nflux; j++ )
1376 {
1377 fref += rfield.flux[0][j]*opac.opacity_abs[j];
1378 fout += rfield.otslin[j]*opac.opacity_abs[j];
1379 fsum += rfield.otscon[j]*opac.opacity_abs[j];
1380 sum += rfield.ConInterOut[j]*opac.opacity_abs[j];
1381 flxin += (rfield.outlin[0][j] + rfield.outlin_noplot[j])*opac.opacity_abs[j];
1382 }
1383 fprintf( save.ipPnunit[ipPun], "%10.2e%10.2e%10.2e%10.2e%10.2e\n",
1384 fref, fout, fsum, sum, flxin );
1385 }
1386 }
1387
1388 else if( strcmp(save.chSave[ipPun],"CONI") == 0 )
1389 {
1390 /* save ionizing continuum */
1391 /* set pointer for possible change in units of energy in continuum
1392 * AnuUnit will give anu in whatever units were set with save units */
1393
1394 if( save.lgSaveEveryZone[ipPun] || (strcmp(chTime,"LAST") == 0) )
1395 {
1396 /* this flag will remember whether we have ever printed anything */
1397 bool lgPrt=false;
1398 if( save.lgSaveEveryZone[ipPun] )
1399 fprintf(save.ipPnunit[ipPun],"#save every zone %li\n", nzone);
1400
1401 /* save ionizing continuum command
1402 * this is option to set lowest energy,
1403 * if no number was entered then this was zero */
1404 if( save.punarg[ipPun][0] <= 0. )
1405 i1 = 1;
1406 else if( save.punarg[ipPun][0] < 100. )
1407 i1 = ipoint(save.punarg[ipPun][0]);
1408 else
1409 i1 = (long int)save.punarg[ipPun][0];
1410
1411 sum = 0.;
1412 for( j=i1-1; j < rfield.nflux; j++ )
1413 {
1414 flxcor = rfield.flux[0][j] +
1415 rfield.otslin[j] +
1416 rfield.otscon[j] +
1417 rfield.ConInterOut[j] +
1418 rfield.outlin[0][j] + rfield.outlin_noplot[j];
1419
1420 sum += flxcor*opac.opacity_abs[j];
1421 }
1422
1423 if( sum > 0. )
1424 sum = 1./sum;
1425 else
1426 sum = 1.;
1427
1428 fsum = 0.;
1429
1430 for( j=i1-1; j<rfield.nflux; ++j)
1431 {
1432 flxcor = rfield.flux[0][j] +
1433 rfield.otslin[j] +
1434 rfield.otscon[j] +
1435 rfield.ConInterOut[j]+
1436 rfield.outlin[0][j] + rfield.outlin_noplot[j];
1437
1438 fsum += flxcor*opac.opacity_abs[j];
1439
1440 /* punched quantities are freq, flux, flux*cross sec,
1441 * fraction of total, integral fraction of total */
1442 RateInter = flxcor*opac.opacity_abs[j]*sum;
1443
1444 /* punage(ipPun,2) is lowest interaction rate to consider, def=0.01 (1 percent) */
1445 /* >>chng 01 nov 22, format to c-friendly */
1446 if( (RateInter >= save.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) )
1447 {
1448 lgPrt = true;
1449 /* >>chng 96 oct 22, format of anu to 11.5 to resolve energy mesh near 1 Ryd */
1450 fprintf( save.ipPnunit[ipPun],
1451 "%li\t%.5e\t%.2e\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2e\t%.2e\t%.4s\t%.4s\n",
1452 j,
1453 AnuUnit(rfield.AnuOrg[j]),
1454 flxcor,
1455 flxcor*opac.opacity_abs[j],
1456 rfield.flux[0][j]/flxcor,
1457 rfield.otslin[j]/flxcor,
1458 rfield.otscon[j]/flxcor,
1459 (rfield.outlin[0][j] + rfield.outlin_noplot[j])/flxcor,
1460 rfield.ConInterOut[j]/flxcor,
1461 RateInter,
1462 fsum*sum,
1463 rfield.chLineLabel[j],
1464 rfield.chContLabel[j] );
1465 }
1466 }
1467 if( !lgPrt )
1468 {
1469 /* entered logical block but did not print anything */
1470 fprintf(save.ipPnunit[ipPun],
1471 " punchdo, the PUNCH IONIZING CONTINUUM command "
1472 "did not find a strong point, sum and fsum were %.2e %.2e\n",
1473 sum,fsum);
1474 fprintf(save.ipPnunit[ipPun],
1475 " punchdo, the low-frequency energy was %.5e Ryd\n",
1476 rfield.anu[i1-1]);
1477 fprintf(save.ipPnunit[ipPun],
1478 " You can reset the threshold for the lowest fractional "
1479 "interaction to print with the second number of the save command\n"
1480 " The fraction was %.3f and this was too large.\n",
1481 save.punarg[ipPun][1]);
1482 }
1483 }
1484 }
1485#ifdef USE_NLTE7
1486 else if( strcmp(save.chSave[ipPun],"CONl") == 0 )
1487 {
1488 if( strcmp(chTime,"LAST") != 0 )
1489 {
1490 int nEmType = (int)save.punarg[ipPun][0];
1491 ASSERT( nEmType==0 || nEmType==1 );
1492
1493 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1494 {
1495
1496 // a value < 0. indicates that energy should be conserved
1497 realnum resolution = ( save.Resolution < realnum(0.) ) ?
1498 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1499
1500 /* the reflected continuum */
1501 flxref = (rfield.anu2[j]*((double)rfield.ConRefIncid[nEmType][j]+rfield.ConEmitReflec[nEmType][j])/rfield.widflx[j] +
1502 rfield.anu[j]*resolution*rfield.reflin[nEmType][j])*EN1RYD;
1503
1504 /* the attenuated incident continuum */
1505 flxatt = rfield.flux[nEmType][j]*rfield.anu2[j]*EN1RYD/
1506 rfield.widflx[j]*radius.r1r0sq *
1507 PrettyTransmission( j, rfield.trans_coef_total[j] );
1508
1509 /* the outward emitted continuum */
1510 conem = (((double)rfield.ConEmitOut[nEmType][j])/
1511 rfield.widflx[j]*rfield.anu2[j] + resolution*
1512 rfield.outlin[nEmType][j]*rfield.anu[j])*radius.r1r0sq*
1513 EN1RYD*geometry.covgeo;
1514
1515 /* sum of emitted and transmitted continua */
1516 flxtrn = conem + flxatt;
1517
1518 //Set upper and lower limits on which wavelength/energy values are printed.
1519 double lowlim, highlim, NRGeV;
1520 lowlim = 1.0;
1521 highlim = 250.0;
1522 NRGeV = AnuUnit(rfield.AnuOrg[j]);
1523
1524 if( NRGeV >= lowlim && NRGeV <= highlim )
1525 {
1526 /* photon energy in appropriate energy or wavelength units */
1527 fprintf( save.ipPnunit[ipPun],"%14.8e ", NRGeV );
1528 /* print zeroes for BB BF and FF */
1529 fprintf( save.ipPnunit[ipPun],"%14.8e %14.8e %14.8e ",0.,0.,0.);
1530 /* total cont */
1531 fprintf( save.ipPnunit[ipPun],"%14.8e\n", (flxref + flxtrn)/NRGeV );
1532 }
1533 }
1534 }
1535
1536 }
1537#endif
1538
1539
1540 else if( strcmp(save.chSave[ipPun],"CONS") == 0 )
1541 {
1542 if( strcmp(chTime,"LAST") != 0 )
1543 {
1544 // continuum volume emissivity and opacity as a function of radius
1545 // command was "save continuum emissivity" */
1546 if( save.ipEmisFreq[ipPun] < 0 )
1547 save.ipEmisFreq[ipPun] = ipoint(save.emisfreq[ipPun].Ryd());
1548 j = save.ipEmisFreq[ipPun]-1;
1549
1550 fprintf( save.ipPnunit[ipPun],
1551 "%.14e\t%.14e\t%.5e\t%.5e\t%.5e\n",
1552 radius.Radius_mid_zone,
1553 radius.depth_mid_zone,
1554 rfield.anu2[j]*rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]*EN1RYD,
1555 opac.opacity_abs[j],
1556 opac.opacity_sct[j] );
1557 }
1558 }
1559
1560 else if( strcmp(save.chSave[ipPun],"CORA") == 0 )
1561 {
1562 /* save raw continuum */
1563 /* set pointer for possible change in units of energy in continuum
1564 * AnuUnit will give anu in whatever units were set with save units */
1565
1566 if( strcmp(chTime,"LAST") == 0 )
1567 {
1568 /* this option to save all raw ionizing continuum */
1569 for( j=0;j<rfield.nflux;j = j + save.ncSaveSkip)
1570 {
1571 fprintf( save.ipPnunit[ipPun],
1572 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t",
1573 AnuUnit(rfield.AnuOrg[j]),
1574 rfield.flux[0][j],
1575 rfield.otslin[j],
1576 rfield.otscon[j],
1577 rfield.ConRefIncid[0][j],
1578 rfield.ConEmitReflec[0][j],
1579 rfield.ConInterOut[j],
1580 rfield.outlin[0][j]+rfield.outlin_noplot[j],
1581 rfield.ConEmitOut[0][j],
1582 rfield.chLineLabel[j],
1583 rfield.chContLabel[j]
1584 );
1585 /* number of lines within that cell */
1586 fprintf( save.ipPnunit[ipPun], "%li\n", rfield.line_count[j] );
1587 }
1588 }
1589 }
1590
1591 else if( strcmp(save.chSave[ipPun],"CONo") == 0 )
1592 {
1593 /* save outward local continuum */
1594 /* set pointer for possible change in units of energy in continuum
1595 * AnuUnit will give anu in whatever units were set with save units */
1596
1597 if( strcmp(chTime,"LAST") == 0 )
1598 {
1599 /* option to save out outward continuum here */
1600 for( j=0;j<rfield.nflux; j = j + save.ncSaveSkip)
1601 {
1602 fprintf( save.ipPnunit[ipPun], "%11.5e%10.2e%10.2e\n",
1603 AnuUnit(rfield.AnuOrg[j]),
1604 rfield.ConEmitOut[0][j]+ rfield.outlin[0][j] + rfield.outlin_noplot[j],
1605 (rfield.flux[0][j] + rfield.otscon[j] + rfield.otslin[j] +
1606 rfield.ConInterOut[j])*opac.opacity_abs[j]*
1607 rfield.anu[j] );
1608 }
1609 }
1610 }
1611
1612 else if( strcmp(save.chSave[ipPun],"CONO") == 0 )
1613 {
1614 /* save outward continuum */
1615 /* set pointer for possible change in units of energy in continuum
1616 * AnuUnit will give anu in whatever units were set with save units */
1617
1618 if( strcmp(chTime,"LAST") == 0 )
1619 {
1620 /* option to save out outward continuum */
1621 for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1622 {
1623 // a value < 0. indicates that energy should be conserved
1624 realnum resolution = ( save.Resolution < realnum(0.) ) ?
1625 rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1626
1627 fprintf( save.ipPnunit[ipPun], "%11.5e%10.2e%10.2e%10.2e%10.2e\n",
1628 AnuUnit(rfield.AnuOrg[j]),
1629 rfield.flux[0][j]*rfield.anu2[j]* EN1RYD/rfield.widflx[j]*radius.r1r0sq,
1630 rfield.ConInterOut[j]/rfield.widflx[j]*rfield.anu2[j]*radius.r1r0sq*EN1RYD,
1631 resolution*(rfield.outlin[0][j]+rfield.outlin_noplot[j])*rfield.anu[j]*radius.r1r0sq*EN1RYD,
1632 (rfield.ConInterOut[j]/rfield.widflx[j]*
1633 rfield.anu2[j] + resolution*(rfield.outlin[0][j]+rfield.outlin_noplot[j])*
1634 rfield.anu[j])*radius.r1r0sq*EN1RYD );
1635 }
1636 }
1637 }
1638
1639 else if( strcmp(save.chSave[ipPun],"CONT") == 0 )
1640 {
1641 /* save transmitted continuum - this is not the main "save continuum"
1642 * command - search on "CON " above
1643 * set pointer for possible change in units of energy in continuum
1644 * AnuUnit will give anu in whatever units were set with save units */
1645
1646 if( strcmp(chTime,"LAST") == 0 )
1647 {
1648 fprintf( save.ipPnunit[ipPun], "#\n" );
1649 fprintf( save.ipPnunit[ipPun], "%32ld # file format version number\n",
1650 VERSION_TRNCON );
1651 fprintf( save.ipPnunit[ipPun], "%s # check 1\n",
1652 continuum.mesh_md5sum.c_str() );
1653 union {
1654 double x;
1655 uint32 i[2];
1656 } u;
1657 u.x = double(rfield.emm);
1658 if( cpu.i().big_endian() )
1659 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 2\n",
1660 u.i[0], u.i[1] );
1661 else
1662 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 2\n",
1663 u.i[1], u.i[0] );
1664 u.x = double(rfield.egamry);
1665 if( cpu.i().big_endian() )
1666 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 3\n",
1667 u.i[0], u.i[1] );
1668 else
1669 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 3\n",
1670 u.i[1], u.i[0] );
1671 u.x = continuum.ResolutionScaleFactor;
1672 if( cpu.i().big_endian() )
1673 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 4\n",
1674 u.i[0], u.i[1] );
1675 else
1676 fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 4\n",
1677 u.i[1], u.i[0] );
1678 fprintf( save.ipPnunit[ipPun], "%32ld # nflux\n",
1679 (rfield.nflux+save.ncSaveSkip-1)/save.ncSaveSkip );
1680 fprintf( save.ipPnunit[ipPun], "#\n" );
1681
1682 const realnum *trans_coef_total=rfield.getCoarseTransCoef();
1683
1684 /* this option to save transmitted continuum */
1685 for( j=0; j < rfield.nflux; j += save.ncSaveSkip )
1686 {
1687 /* attenuated incident continuum
1688 * >>chng 97 jul 10, remove SaveLWidth from this one only since
1689 * we must conserve energy even in lines
1690 * >>chng 07 apr 26 include transmission coefficient */
1691 flxatt = rfield.flux[0][j]*rfield.anu2[j]*EN1RYD/
1692 rfield.widflx[j]*radius.r1r0sq*trans_coef_total[j];
1693
1694 /*conem = (rfield.ConOutNoInter[j] + rfield.ConInterOut[j]+rfield.outlin[0][j])*
1695 rfield.anu2[j];
1696 conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;*/
1697 /* >>chng 00 jan 03, above did not include all contributors.
1698 * Pasted in below from usual
1699 * save continuum command */
1700 /* >>chng 04 jul 15, removed factor of save.SaveLWidth -
1701 * this should not be there to conserve energy, as explained in hazy
1702 * where command was documented, and in comment above. caught by PvH */
1703 /* >>chng 04 jul 23, incorrect use of outlin - before multiplied by an2,
1704 * quantity should be photons per Ryd, since init quantity is
1705 * photons per cell. Must div by widflx. caught by PvH */
1706 /*conem = (rfield.ConEmitOut[0][j]/rfield.widflx[j]*rfield.anu2[j] +
1707 rfield.outlin[0][j]*rfield.anu[j])*radius.r1r0sq*
1708 EN1RYD*geometry.covgeo;*/
1709 conem = (rfield.ConEmitOut[0][j] + rfield.outlin[0][j]) / rfield.widflx[j]*
1710 rfield.anu2[j]*radius.r1r0sq*EN1RYD*geometry.covgeo;
1711
1712 flxtrn = conem + flxatt;
1713
1714 /* use AnuOrg here instead of anu since probably
1715 * going to be used by table read
1716 * and we want good anu array for sanity check*/
1717 fprintf( save.ipPnunit[ipPun],"%.5e\t%.3e\t%.3e\n",
1718 AnuUnit(rfield.AnuOrg[j]), flxtrn,
1719 trans_coef_total[j] );
1720 }
1721 }
1722 }
1723
1724 else if( strcmp(save.chSave[ipPun],"CON2") == 0 )
1725 {
1726 /* save total two-pohoton continuum */
1727 if( strcmp(chTime,"LAST") == 0 )
1728 {
1729 /* this option to save diffuse continuum */
1730 for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1731 {
1732 fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\n",
1733 AnuUnit(rfield.AnuOrg[j]),
1734 rfield.TotDiff2Pht[j]/rfield.widflx[j] ,
1735 rfield.TotDiff2Pht[j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
1736 }
1737 }
1738 }
1739
1740 else if( strcmp(save.chSave[ipPun],"DUSE") == 0 )
1741 {
1742 /* save grain extinction - includes only grain opacity, not total */
1743 if( strcmp(chTime,"LAST") != 0 )
1744 {
1745 fprintf( save.ipPnunit[ipPun], " %.5e\t",
1746 radius.depth_mid_zone );
1747
1748 /* visual extinction of an extended source (like a PDR)*/
1749 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
1750
1751 /* visual extinction of point source (star)*/
1752 fprintf( save.ipPnunit[ipPun], "%.2e\n" , rfield.extin_mag_V_point);
1753 }
1754 }
1755
1756 else if( strcmp(save.chSave[ipPun],"DUSO") == 0 )
1757 {
1758 /* save grain opacity */
1759 if( strcmp(chTime,"LAST") == 0 )
1760 {
1761 for( j=0; j < rfield.nflux; j++ )
1762 {
1763 double scat;
1764 fprintf( save.ipPnunit[ipPun],
1765 "%.5e\t%.2e\t%.2e\t%.2e\t",
1766 /* photon energy or wavelength */
1767 AnuUnit(rfield.AnuOrg[j]),
1768 /* total opacity, discount forward scattering */
1769 gv.dstab[j] + gv.dstsc[j],
1770 /* absorption opacity */
1771 gv.dstab[j],
1772 /* scatter, with forward discounted */
1773 gv.dstsc[j] );
1774 /* add together total scattering, discounting 1-g */
1775 scat = 0.;
1776 /* sum over all grain species */
1777 for( size_t nd=0; nd < gv.bin.size(); nd++ )
1778 {
1779 scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund;
1780 }
1781 /* finally, scattering including effects of forward scattering */
1782 fprintf( save.ipPnunit[ipPun],
1783 "%.2e\t", scat );
1784 fprintf( save.ipPnunit[ipPun],
1785 "%.2e\n", gv.dstsc[j]/SDIV(gv.dstab[j] + gv.dstsc[j]) );
1786 }
1787 }
1788 }
1789
1790 /* save grain abundance and save grain D/G ratio commands */
1791 else if( strcmp(save.chSave[ipPun],"DUSA") == 0 ||
1792 strcmp(save.chSave[ipPun],"DUSD") == 0 )
1793 {
1794 bool lgDGRatio = ( strcmp(save.chSave[ipPun],"DUSD") == 0 );
1795
1796 /* grain abundance */
1797 if( strcmp(chTime,"LAST") != 0 )
1798 {
1799 /* used to print header exactly one time */
1800 static bool lgMustPrtHeaderDRRatio = true,
1801 lgMustPrtHeaderAbundance=true;
1802 /* print grain header first if this has not yet been done */
1803 if( ( lgMustPrtHeaderDRRatio && lgDGRatio ) ||
1804 ( lgMustPrtHeaderAbundance && !lgDGRatio ) )
1805 {
1806 /* only print one header for each case, but must
1807 * track separately if both used in same sim */
1808 if( lgMustPrtHeaderDRRatio && lgDGRatio )
1809 lgMustPrtHeaderDRRatio = false;
1810 else if( lgMustPrtHeaderAbundance &&!lgDGRatio )
1811 lgMustPrtHeaderAbundance = false;
1812
1813 fprintf( save.ipPnunit[ipPun], "#Depth" );
1814 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1815 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1816 fprintf( save.ipPnunit[ipPun], "\ttotal\n" );
1817
1818 /* now print grain radius converting from cm to microns */
1819 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1820 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1821 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1822 fprintf( save.ipPnunit[ipPun], "\txxxx\n" );
1823 }
1824 fprintf( save.ipPnunit[ipPun], " %.5e",
1825 radius.depth_mid_zone );
1826 /* grain abundance per bin in g/cm^3 */
1827 double total = 0.;
1828 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1829 {
1830 double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]*
1831 gv.bin[nd]->cnv_H_pCM3;
1832 if( lgDGRatio )
1833 abund /= dense.xMassDensity;
1834 fprintf( save.ipPnunit[ipPun], "\t%.3e", abund );
1835 total += abund;
1836 }
1837 fprintf( save.ipPnunit[ipPun], "\t%.3e\n", total );
1838 }
1839 }
1840
1841 else if( strcmp(save.chSave[ipPun],"DUSP") == 0 )
1842 {
1843 /* grain potential */
1844 if( strcmp(chTime,"LAST") != 0 )
1845 {
1846 /* used to print header exactly one time */
1847 static bool lgMustPrtHeader = true;
1848 /* do labels first if this is first zone */
1849 if( lgMustPrtHeader )
1850 {
1851 /* first print string giving grain id */
1852 fprintf( save.ipPnunit[ipPun], "#Depth" );
1853 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1854 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1855 fprintf( save.ipPnunit[ipPun], "\n" );
1856
1857 /* now print grain radius converting from cm to microns */
1858 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1859 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1860 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1861 fprintf( save.ipPnunit[ipPun], "\n" );
1862
1863 /* don't need to do this, ever again */
1864 lgMustPrtHeader = false;
1865 }
1866 fprintf( save.ipPnunit[ipPun], " %.5e",
1867 radius.depth_mid_zone );
1868 /* grain potential in eV */
1869 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1870 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD );
1871 fprintf( save.ipPnunit[ipPun], "\n" );
1872 }
1873 }
1874
1875 else if( strcmp(save.chSave[ipPun],"DUSR") == 0 )
1876 {
1877 /* grain H2 formation rates */
1878 if( strcmp(chTime,"LAST") != 0 )
1879 {
1880 /* used to print header exactly one time */
1881 static bool lgMustPrtHeader = true;
1882 /* do labels first if this is first zone */
1883 if( lgMustPrtHeader )
1884 {
1885 /* first print string giving grain id */
1886 fprintf( save.ipPnunit[ipPun], "#Depth" );
1887 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1888 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1889 fprintf( save.ipPnunit[ipPun], "\n" );
1890
1891 /* now print grain radius converting from cm to microns */
1892 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1893 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1894 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1895 fprintf( save.ipPnunit[ipPun], "\n" );
1896
1897 /* don't need to do this, ever again */
1898 lgMustPrtHeader = false;
1899 }
1900 fprintf( save.ipPnunit[ipPun], " %.5e",
1901 radius.depth_mid_zone );
1902 /* grain formation rate for H2 */
1903 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1904 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used );
1905 fprintf( save.ipPnunit[ipPun], "\n" );
1906 }
1907 }
1908
1909 else if( strcmp(save.chSave[ipPun],"DUST") == 0 )
1910 {
1911 /* grain temperatures - K*/
1912 if( strcmp(chTime,"LAST") != 0 )
1913 {
1914 /* used to print header exactly one time */
1915 static bool lgMustPrtHeader = true;
1916 /* do labels first if this is first zone */
1917 if( lgMustPrtHeader )
1918 {
1919 /* first print string giving grain id */
1920 fprintf( save.ipPnunit[ipPun], "#Depth" );
1921 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1922 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1923 fprintf( save.ipPnunit[ipPun], "\n" );
1924
1925 /* now print grain radius converting from cm to microns */
1926 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1927 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1928 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1929 fprintf( save.ipPnunit[ipPun], "\n" );
1930
1931 /* don't need to do this, ever again */
1932 lgMustPrtHeader = false;
1933 }
1934 fprintf( save.ipPnunit[ipPun], " %.5e",
1935 radius.depth_mid_zone );
1936 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1937 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust );
1938 fprintf( save.ipPnunit[ipPun], "\n" );
1939 }
1940 }
1941
1942 else if( strcmp(save.chSave[ipPun],"DUSC") == 0 )
1943 {
1944 /* save grain charge - eden from grains and
1945 * charge per grain in electrons / grain */
1946 if( strcmp(chTime,"LAST") != 0 )
1947 {
1948 /* used to print header exactly one time */
1949 static bool lgMustPrtHeader = true;
1950 /* do labels first if this is first zone */
1951 if( lgMustPrtHeader )
1952 {
1953 /* first print string giving grain id */
1954 fprintf( save.ipPnunit[ipPun], "#Depth\tne(grn)" );
1955 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1956 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1957 fprintf( save.ipPnunit[ipPun], "\n" );
1958
1959 /* now print grain radius converting from cm to microns */
1960 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)\tne\t" );
1961 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1962 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1963 fprintf( save.ipPnunit[ipPun], "\n" );
1964
1965 /* don't need to do this, ever again */
1966 lgMustPrtHeader = false;
1967 }
1968
1969 fprintf( save.ipPnunit[ipPun], " %.5e\t%.4e",
1970 radius.depth_mid_zone ,
1971 /* electron density contributed by grains, in e/cm^3,
1972 * positive number means grain supplied free electrons */
1973 gv.TotalEden );
1974
1975 /* average charge per grain in electrons */
1976 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1977 {
1978 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ );
1979 }
1980 fprintf( save.ipPnunit[ipPun], "\n" );
1981 }
1982 }
1983
1984 else if( strcmp(save.chSave[ipPun],"DUSH") == 0 )
1985 {
1986 /* grain heating */
1987 if( strcmp(chTime,"LAST") != 0 )
1988 {
1989 /* used to print header exactly one time */
1990 static bool lgMustPrtHeader = true;
1991 /* save grain charge, but do labels first if this is first zone */
1992 if( lgMustPrtHeader )
1993 {
1994 /* first print string giving grain id */
1995 fprintf( save.ipPnunit[ipPun], "#Depth" );
1996 for( size_t nd=0; nd < gv.bin.size(); ++nd )
1997 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1998 fprintf( save.ipPnunit[ipPun], "\n" );
1999
2000 /* now print grain radius converting from cm to microns */
2001 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
2002 for( size_t nd=0; nd < gv.bin.size(); ++nd )
2003 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
2004 fprintf( save.ipPnunit[ipPun], "\n" );
2005
2006 /* don't need to do this, ever again */
2007 lgMustPrtHeader = false;
2008 }
2009 fprintf( save.ipPnunit[ipPun], " %.5e",
2010 radius.depth_mid_zone );
2011 /* grain heating */
2012 for( size_t nd=0; nd < gv.bin.size(); ++nd )
2013 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl );
2014 fprintf( save.ipPnunit[ipPun], "\n" );
2015 }
2016 }
2017
2018 else if( strcmp(save.chSave[ipPun],"DUSV") == 0 )
2019 {
2020 /* grain drift velocities */
2021 if( strcmp(chTime,"LAST") != 0 )
2022 {
2023 /* used to print header exactly one time */
2024 static bool lgMustPrtHeader = true;
2025 /* save grain velocity, but do labels first if this is first zone */
2026 if( lgMustPrtHeader )
2027 {
2028 /* first print string giving grain id */
2029 fprintf( save.ipPnunit[ipPun], "#Depth" );
2030 for( size_t nd=0; nd < gv.bin.size(); ++nd )
2031 fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
2032 fprintf( save.ipPnunit[ipPun], "\n" );
2033
2034 /* now print grain radius converting from cm to microns */
2035 fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
2036 for( size_t nd=0; nd < gv.bin.size(); ++nd )
2037 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
2038 fprintf( save.ipPnunit[ipPun], "\n" );
2039
2040 /* don't need to do this, ever again */
2041 lgMustPrtHeader = false;
2042 }
2043 fprintf( save.ipPnunit[ipPun], " %.5e",
2044 radius.depth_mid_zone );
2045 /* grain drift velocity in km/s */
2046 for( size_t nd=0; nd < gv.bin.size(); ++nd )
2047 fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 );
2048 fprintf( save.ipPnunit[ipPun], "\n" );
2049 }
2050 }
2051
2052 /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
2053 else if( strcmp(save.chSave[ipPun],"DUSQ") == 0 )
2054 {
2055 /* save grain Qs */
2056 if( strcmp(chTime,"LAST") == 0 )
2057 {
2058 for( j=0; j < rfield.nflux; j++ )
2059 {
2060 fprintf( save.ipPnunit[ipPun], "%11.4e",
2061 rfield.anu[j] );
2062 for( size_t nd=0; nd < gv.bin.size(); nd++ )
2063 {
2064 fprintf( save.ipPnunit[ipPun], "%9.1e%9.1e",
2065 gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea,
2066 gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea );
2067 }
2068 fprintf( save.ipPnunit[ipPun], "\n" );
2069 }
2070 }
2071 }
2072
2073 else if( strcmp(save.chSave[ipPun],"ELEM") == 0 )
2074 {
2075 if( strcmp(chTime,"LAST") != 0 )
2076 {
2077 realnum renorm = 1.f;
2078
2079 /* this is the index for the atomic number on the physical scale */
2080 /* >>chng 04 nov 23, use c scale throughout */
2081 nelem = (long int)save.punarg[ipPun][0];
2082 ASSERT( nelem >= ipHYDROGEN );
2083
2084 /* don't do this if element is not turned on */
2085 if( dense.lgElmtOn[nelem] )
2086 {
2087 /* >>chng 04 nov 23, add density option, leave as cm-3
2088 * default is still norm to total of that element */
2089 if( save.punarg[ipPun][1] == 0 )
2090 renorm = dense.gas_phase[nelem];
2091
2092 fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone );
2093
2094 for( j=0; j <= (nelem + 1); ++j)
2095 {
2096 fprintf( save.ipPnunit[ipPun], "\t%.2e",
2097 dense.xIonDense[nelem][j]/renorm );
2098 }
2099 if( nelem==ipHYDROGEN )
2100 {
2101 /* H2 */
2102 fprintf( save.ipPnunit[ipPun], "\t%.2e",
2103 hmi.H2_total/renorm );
2104 }
2105 /* >>chng 04 nov 23 add C and O fine structure pops */
2106 else if( nelem==ipCARBON )
2107 {
2108 fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
2109 colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm);
2110 fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e",
2111 colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm);
2112 fprintf( save.ipPnunit[ipPun], "\t%.2e",
2113 findspecieslocal("CO")->den/renorm );
2114 }
2115 else if( nelem==ipOXYGEN )
2116 {
2117 fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
2118 colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm);
2119 }
2120 fprintf( save.ipPnunit[ipPun], "\n" );
2121 }
2122 }
2123 }
2124
2125 else if( strcmp(save.chSave[ipPun],"RECA") == 0 )
2126 {
2127 /* this will create table for AGN3 then exit,
2128 * routine is in makerecom.c */
2129 ion_recombAGN( save.ipPnunit[ipPun] );
2130 cdEXIT(EXIT_FAILURE);
2131 }
2132
2133 else if( strcmp(save.chSave[ipPun],"RECE") == 0 )
2134 {
2135 /* save recombination efficiencies,
2136 * option turned on with the "save recombination efficiencies" command
2137 * output for the save recombination coefficients command is actually
2138 * produced by a series of routines, as they generate the recombination
2139 * coefficients. these include
2140 * dielsupres, helium, hydrorecom, iibod, and makerecomb*/
2141 fprintf( save.ipPnunit[ipPun],
2142 "%12.4e %12.4e %12.4e %12.4e\n",
2143 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecRad],
2144 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecNetEsc] ,
2145 iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecRad],
2146 iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecNetEsc]);
2147 }
2148
2149 else if( strncmp(save.chSave[ipPun],"FEc" , 3 ) == 0 )
2150 {
2151 /* FeII continuum */
2152 if( strcmp(chTime,"LAST") == 0 )
2153 {
2154
2155 if( save.punarg[ipPun][0] == 1 )
2156 {
2157 // row format used for grids - one time print of wavelength grid first
2158 char chTempHeader[100];
2159 long ipFeII_Cont_type;
2160 if( strcmp(save.chSave[ipPun],"FEcI") == 0 )
2161 {
2162 strcpy(chTempHeader , "#FeII inward: Wl(A)\tInt[erg cm-2 s-1]\n");
2163 ipFeII_Cont_type = 1;
2164 }
2165 else if( strcmp(save.chSave[ipPun],"FEcO") == 0 )
2166 {
2167 strcpy(chTempHeader , "#FeII outward: Wl(A)\tInt[erg cm-2 s-1]\n");
2168 ipFeII_Cont_type = 2;
2169 }
2170 else if( strcmp(save.chSave[ipPun],"FEcT") == 0 )
2171 {
2172 strcpy(chTempHeader , "#FeII total: Wl(A)\tInt[erg cm-2 s-1]\n");
2173 ipFeII_Cont_type = 3;
2174 }
2175 else
2176 TotalInsanity();
2177
2178 if( save.lgPunHeader[ipPun] &&
2179 (!grid.lgGrid || optimize.nOptimiz==0) )
2180 {
2181 // one time print of FeII continuum header
2182 j = 0;
2183 fprintf( save.ipPnunit[ipPun], "%s%.5f",
2184 chTempHeader ,
2185 AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
2186 for( j=1; j < nFeIIConBins; j++ )
2187 {
2188 fprintf( save.ipPnunit[ipPun], "\t%.5e",
2189 AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
2190 }
2191 fprintf( save.ipPnunit[ipPun], "\n");
2192 save.lgPunHeader[ipPun] = false;
2193 }
2194 // give intensities
2195 fprintf( save.ipPnunit[ipPun], "%.2f",
2196 SaveFeII_cont( 0 , ipFeII_Cont_type ));
2197 for( j=1; j < nFeIIConBins; j++ )
2198 {
2199 fprintf( save.ipPnunit[ipPun], "\t%e",
2200 SaveFeII_cont( j , ipFeII_Cont_type ) );
2201 }
2202 fprintf( save.ipPnunit[ipPun], "\n");
2203 }
2204 else if( save.punarg[ipPun][0] == 2 )
2205 {
2206 // the default, four columns, wl, total, inward, outward
2207 // one time print of header
2208 if( save.lgPunHeader[ipPun] &&
2209 (!grid.lgGrid || optimize.nOptimiz==0) )
2210 {
2211 fprintf( save.ipPnunit[ipPun], "FeiI wl, total, inward, outward\n");
2212 save.lgPunHeader[ipPun] = false;
2213 }
2214 for( j=0; j < nFeIIConBins; j++ )
2215 {
2216 fprintf( save.ipPnunit[ipPun], "%.5f",
2217 AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
2218 fprintf( save.ipPnunit[ipPun], "\t%e",
2219 SaveFeII_cont( j , 3 ) );// total
2220 fprintf( save.ipPnunit[ipPun], "\t%e",
2221 SaveFeII_cont( j , 1 ) );// inward
2222 fprintf( save.ipPnunit[ipPun], "\t%e",
2223 SaveFeII_cont( j , 2 ) );// outward
2224 fprintf( save.ipPnunit[ipPun], "\n");
2225 }
2226 }
2227 }
2228 }
2229
2230 /* save column densities */
2231 else if( strcmp(save.chSave[ipPun],"FENl") == 0 )
2232 {
2233 if( strcmp(chTime,"LAST") == 0 )
2234 {
2235 /* save FeII level energies and stat weights, followed by column density */
2236 FeIIPunchColden( save.ipPnunit[ipPun] );
2237 }
2238 }
2239
2240 else if( strcmp(save.chSave[ipPun],"FE2l") == 0 )
2241 {
2242 if( strcmp(chTime,"LAST") == 0 )
2243 {
2244 /* save FeII level energies and stat weights */
2245 FeIIPunchLevels( save.ipPnunit[ipPun] );
2246 }
2247 }
2248
2249 else if( strcmp(save.chSave[ipPun],"FEli") == 0 )
2250 {
2251 if( strcmp(chTime,"LAST") == 0 )
2252 {
2253 /* save line intensities, routine is in atom_feii.c */
2254 FeIISaveLines( save.ipPnunit[ipPun] );
2255 }
2256 }
2257
2258 else if( strcmp(save.chSave[ipPun],"FE2o") == 0 )
2259 {
2260 if( strcmp(chTime,"LAST") == 0 )
2261 {
2262 /* save line optical depths, routine is in atom_feii.c */
2263 FeIIPunchOpticalDepth( save.ipPnunit[ipPun] );
2264 }
2265 }
2266
2267 else if( strcmp(save.chSave[ipPun],"FRED") == 0 )
2268 {
2269 /* set with save Fred command, this punches some stuff from
2270 * Fred Hamann's dynamics project */
2271 if( strcmp(chTime,"LAST") != 0 )
2272 {
2273 /* Fred's list */
2274 fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.3e\t%.3e\t%.3e"
2275 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2276 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2277 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2278 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2279 "\t%.3e\t%.3e\n",
2280 radius.Radius, radius.depth ,wind.windv/1e5,
2281 wind.dvdr,
2282 dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te,
2283 wind.AccelLine , wind.AccelCont ,
2284 wind.fmul ,
2285 // acceleration in this zone due to electron scattering,
2286 // if incident SED was not attenuated
2287 pressure.pinzon_PresIntegElecThin/dense.xMassDensity/radius.drad_x_fillfac ,
2288 mean.xIonMean[0][ipHYDROGEN][0][0] , mean.xIonMean[0][ipHYDROGEN][1][0] ,
2289 mean.xIonMean[0][ipHELIUM][0][0] , mean.xIonMean[0][ipHELIUM][1][0] ,
2290 mean.xIonMean[0][ipHELIUM][2][0] ,
2291 mean.xIonMean[0][ipCARBON][1][0] , mean.xIonMean[0][ipCARBON][2][0] ,
2292 mean.xIonMean[0][ipCARBON][3][0] ,
2293 mean.xIonMean[0][ipOXYGEN][0][0] , mean.xIonMean[0][ipOXYGEN][1][0] ,
2294 mean.xIonMean[0][ipOXYGEN][2][0] , mean.xIonMean[0][ipOXYGEN][3][0] ,
2295 mean.xIonMean[0][ipOXYGEN][4][0] , mean.xIonMean[0][ipOXYGEN][5][0] ,
2296 mean.xIonMean[0][ipOXYGEN][6][0] , mean.xIonMean[0][ipOXYGEN][7][0] ,
2297 dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] ,
2298 dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] ,
2299 dense.xIonDense[ipHELIUM][2] ,
2300 dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] ,
2301 dense.xIonDense[ipCARBON][3] ,
2302 dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] ,
2303 dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] ,
2304 dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] ,
2305 dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] ,
2306 mean.xIonMean[0][ipMAGNESIUM][1][0] , dense.xIonDense[ipMAGNESIUM][1],
2307 TauLines[ipT1032].Emis().TauIn() ,
2308 TauLines[ipT1032].Emis().TauCon() );
2309 }
2310 }
2311
2312 else if( strcmp(save.chSave[ipPun],"FE2d") == 0 )
2313 {
2314 /* save some departure coefficients for large FeII atom */
2315 if( strcmp(chTime,"LAST") != 0 )
2316 FeIIPunDepart(save.ipPnunit[ipPun],false);
2317 }
2318
2319 else if( strcmp(save.chSave[ipPun],"FE2D") == 0 )
2320 {
2321 /* save all departure coefficients for large FeII atom */
2322 if( strcmp(chTime,"LAST") != 0 )
2323 FeIIPunDepart(save.ipPnunit[ipPun],true);
2324 }
2325
2326 else if( strcmp(save.chSave[ipPun],"FE2p") == 0 )
2327 {
2328 bool lgFlag = false;
2329 if( save.punarg[ipPun][2] )
2330 lgFlag = true;
2331 /* save small subset of level populations for large FeII atom */
2332 if( strcmp(chTime,"LAST") != 0 )
2333 FeIIPunPop(save.ipPnunit[ipPun],false,0,0,
2334 lgFlag);
2335 }
2336
2337 else if( strcmp(save.chSave[ipPun],"FE2P") == 0 )
2338 {
2339 bool lgFlag = false;
2340 if( save.punarg[ipPun][2] )
2341 lgFlag = true;
2342 /* save range of level populations for large FeII atom */
2343 if( strcmp(chTime,"LAST") != 0 )
2344 FeIIPunPop(save.ipPnunit[ipPun],
2345 true,
2346 (long int)save.punarg[ipPun][0],
2347 (long int)save.punarg[ipPun][1],
2348 lgFlag);
2349 }
2350
2351 /* save spectra in fits format */
2352 else if( strcmp(save.chSave[ipPun],"FITS") == 0 )
2353 {
2354 if( strcmp(chTime,"LAST") == 0 )
2355 {
2356 saveFITSfile( save.ipPnunit[ipPun], NUM_OUTPUT_TYPES );
2357 }
2358 }
2359 /* save gammas (but without element) */
2360 else if( strcmp(save.chSave[ipPun],"GAMt") == 0 )
2361 {
2362 if( strcmp(chTime,"LAST") != 0 )
2363 {
2364 long ns;
2365 /* save photoionization rates, with the PUNCH GAMMAS command */
2366 for( nelem=0; nelem < LIMELM; nelem++ )
2367 {
2368 for( ion=0; ion <= nelem; ion++ )
2369 {
2370 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
2371 {
2372 fprintf( save.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e",
2373 nelem+1, ion+1, ns+1,
2374 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
2375 ionbal.PhotoRate_Shell[nelem][ion][ns][1] ,
2376 ionbal.PhotoRate_Shell[nelem][ion][ns][2] );
2377
2378 for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ )
2379 {
2380 fprintf( save.ipPnunit[ipPun], "%5.2f",
2381 t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) );
2382 }
2383 fprintf( save.ipPnunit[ipPun], "\n" );
2384 }
2385 }
2386 }
2387 }
2388 }
2389
2390 /* save gammas element, ion */
2391 else if( strcmp(save.chSave[ipPun],"GAMe") == 0 )
2392 {
2393 if( strcmp(chTime,"LAST") != 0 )
2394 {
2395 int ns;
2396 nelem = (long)save.punarg[ipPun][0];
2397 ion = (long)save.punarg[ipPun][1];
2398 /* valence shell */
2399 ns = Heavy.nsShells[nelem][ion]-1;
2400 /* show what some of the ionization sources are */
2401 GammaPrt(
2402 opac.ipElement[nelem][ion][ns][0] ,
2403 opac.ipElement[nelem][ion][ns][1] ,
2404 opac.ipElement[nelem][ion][ns][2] ,
2405 save.ipPnunit[ipPun],
2406 ionbal.PhotoRate_Shell[nelem][ion][ns][0] ,
2407 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 );
2408 }
2409 }
2410
2411 else if( strcmp(save.chSave[ipPun],"GAUN") == 0 )
2412 {
2413 /* save gaunt factors */
2414 if( strcmp(chTime,"LAST") != 0 )
2415 SaveGaunts(save.ipPnunit[ipPun]);
2416 }
2417
2418 // generating the SAVE GRID output has been moved to cdPrepareExit()
2419 // to make sure that the output always records any type of failure
2420 }
2421 else
2422 {
2423 //no hit this branch, key should be in next
2424 lgNoHitFirstBranch = true;
2425 }
2426
2427 // hack needed for code to compile with Visual Studio
2428 // keep this identical to the if-statement further up!!
2429 if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
2430 ( lgAbort && strcmp(chTime,"LAST") == 0 ) ||
2431 ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) )
2432 {
2433 if( strcmp(save.chSave[ipPun],"HISp") == 0 )
2434 {
2435 /* save pressure history of current zone */
2436 if( strcmp(chTime,"LAST") != 0 )
2437 {
2438 /* note if pressure convergence failure occurred in history that follows */
2439 if( !conv.lgConvPres )
2440 {
2441 fprintf( save.ipPnunit[ipPun],
2442 "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n",
2443 iteration , nzone );
2444 }
2445 /* note if temperature convergence failure occurred in history that follows */
2446 if( !conv.lgConvTemp )
2447 {
2448 fprintf( save.ipPnunit[ipPun],
2449 "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n",
2450 iteration , nzone );
2451 }
2452 for( unsigned long k=0; k < conv.hist_pres_density.size(); ++k )
2453 {
2454 /* save history of density - pressure, with correct pressure */
2455 fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
2456 iteration,
2457 nzone,
2458 conv.hist_pres_density[k],
2459 conv.hist_pres_current[k],
2460 conv.hist_pres_error[k]);
2461 }
2462 }
2463 }
2464
2465 else if( strcmp(save.chSave[ipPun],"HISt") == 0 )
2466 {
2467 /* save temperature history of current zone */
2468 if( strcmp(chTime,"LAST") != 0 )
2469 {
2470 /* note if pressure convergence failure occurred in history that follows */
2471 if( !conv.lgConvPres )
2472 {
2473 fprintf( save.ipPnunit[ipPun],
2474 "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n",
2475 iteration , nzone );
2476 }
2477 /* note if temperature convergence failure occurred in history that follows */
2478 if( !conv.lgConvTemp )
2479 {
2480 fprintf( save.ipPnunit[ipPun],
2481 "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n",
2482 iteration , nzone );
2483 }
2484 for( unsigned long k=0; k < conv.hist_temp_temp.size(); ++k )
2485 {
2486 /* save history of density - pressure, with correct pressure */
2487 fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
2488 iteration,
2489 nzone,
2490 conv.hist_temp_temp[k],
2491 conv.hist_temp_heat[k],
2492 conv.hist_temp_cool[k]);
2493 }
2494 }
2495 }
2496
2497 else if( strncmp(save.chSave[ipPun],"H2",2) == 0 )
2498 {
2499 /* all save info on large H2 molecule include H2 PDR pdr */
2500 save.whichDiatomToPrint[ipPun]->H2_PunchDo( save.ipPnunit[ipPun] , save.chSave[ipPun] , chTime, ipPun );
2501 }
2502
2503 else if( strcmp(save.chSave[ipPun],"HEAT") == 0 )
2504 {
2505 /* save heating, routine in file of same name */
2506 if( strcmp(chTime,"LAST") != 0 )
2507 SaveHeat(save.ipPnunit[ipPun]);
2508 }
2509
2510 else if( strncmp(save.chSave[ipPun],"HE",2) == 0 )
2511 {
2512 /* various save helium commands */
2513 /* save helium line wavelengths */
2514 if( strcmp(save.chSave[ipPun] , "HELW") == 0 )
2515 {
2516 if( strcmp(chTime,"LAST") == 0 )
2517 {
2518 /* save helium & he-like wavelengths, first header */
2519 fprintf( save.ipPnunit[ipPun],
2520 "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S"
2521 "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" );
2522 fprintf( save.ipPnunit[ipPun], "\n" );
2523 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
2524 {
2525 /* print element name, nuclear charge */
2526 fprintf( save.ipPnunit[ipPun], "%li\t%s",
2527 nelem+1 , elementnames.chElementSym[nelem] );
2528 /*prt_wl print floating wavelength in Angstroms, in output format */
2529 fprintf( save.ipPnunit[ipPun], "\t" );
2530 prt_wl( save.ipPnunit[ipPun] ,
2531 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).WLAng() );
2532 fprintf( save.ipPnunit[ipPun], "\t" );
2533 prt_wl( save.ipPnunit[ipPun] ,
2534 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe1s1S).WLAng() );
2535 fprintf( save.ipPnunit[ipPun], "\t" );
2536 prt_wl( save.ipPnunit[ipPun] ,
2537 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe1s1S).WLAng() );
2538 fprintf( save.ipPnunit[ipPun], "\t" );
2539 prt_wl( save.ipPnunit[ipPun] ,
2540 iso_sp[ipHE_LIKE][nelem].trans(ipHe2s3S,ipHe1s1S).WLAng() );
2541 fprintf( save.ipPnunit[ipPun], "\t" );
2542 prt_wl( save.ipPnunit[ipPun] ,
2543 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe2s3S).WLAng() );
2544 fprintf( save.ipPnunit[ipPun], "\t" );
2545 prt_wl( save.ipPnunit[ipPun] ,
2546 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe2s3S).WLAng() );
2547 fprintf( save.ipPnunit[ipPun], "\t" );
2548 prt_wl( save.ipPnunit[ipPun] ,
2549 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P0,ipHe2s3S).WLAng() );
2550 fprintf( save.ipPnunit[ipPun], "\n");
2551 }
2552 }
2553 }
2554 else
2555 TotalInsanity();
2556 }
2557
2558 /* save hummer, results needed for Lya transport, to feed into David's routine */
2559 else if( strcmp(save.chSave[ipPun],"HUMM") == 0 )
2560 {
2561 eps = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul()/
2562 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Coll().ColUL( colliders );
2563 fprintf( save.ipPnunit[ipPun],
2564 " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
2565 radius.depth_mid_zone,
2566 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2567 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(),
2568 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop(),
2569 phycon.te, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().damp(), eps );
2570 }
2571
2572 else if( strncmp( save.chSave[ipPun] , "HYD", 3 ) == 0 )
2573 {
2574 /* various save hydrogen commands */
2575 if( strcmp(save.chSave[ipPun],"HYDc") == 0 )
2576 {
2577 if( strcmp(chTime,"LAST") != 0 )
2578 {
2579 /* save hydrogen physical conditions */
2580 fprintf( save.ipPnunit[ipPun],
2581 " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2582 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden,
2583 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
2584 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
2585 hmi.H2_total/dense.gas_phase[ipHYDROGEN],
2586 findspecieslocal("H2+")->den/dense.gas_phase[ipHYDROGEN],
2587 findspecieslocal("H3+")->den/dense.gas_phase[ipHYDROGEN],
2588 findspecieslocal("H-")->den/dense.gas_phase[ipHYDROGEN] );
2589 }
2590 }
2591
2592 else if( strcmp(save.chSave[ipPun],"HYDi") == 0 )
2593 {
2594 if( strcmp(chTime,"LAST") != 0 )
2595 {
2596 /* save hydrogen ionization
2597 * this will be total decays to ground */
2598 RateInter = 0.;
2599 stage = iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ColIoniz*dense.eden*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2600 fref = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2601 fout = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2602 /* 06 aug 28, from numLevels_max to _local. */
2603 for( ion=ipH2s; ion < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local; ion++ )
2604 {
2605 /* this is total decays to ground */
2606 RateInter +=
2607 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Aul()*
2608 (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pesc() +
2609 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pelec_esc() +
2610 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pdest());
2611 /* total photo from all levels */
2612 fref += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2613 /* total col ion from all levels */
2614 stage += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].ColIoniz*dense.eden*
2615 iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2616 fout += iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2617 }
2618
2619 /* make these relative to parent ion */
2620 stage /= dense.xIonDense[ipHYDROGEN][1];
2621 fref /= dense.xIonDense[ipHYDROGEN][1];
2622 fout /= dense.xIonDense[ipHYDROGEN][1];
2623
2624 fprintf( save.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e",
2625 nzone,
2626 /* photo and collision ion rates have units s-1 */
2627 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,
2628 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ColIoniz* dense.EdenHCorr,
2629 ionbal.RateRecomTot[ipHYDROGEN][0] );
2630
2631 fprintf( save.ipPnunit[ipPun], "\t%.2e",
2632 iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB );
2633
2634 fprintf( save.ipPnunit[ipPun],
2635 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2636 dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0],
2637 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc/(ionbal.RateRecomTot[ipHYDROGEN][0]),
2638 iso_sp[ipH_LIKE][ipHYDROGEN].fb[1].RadRecomb[ipRecEsc],
2639 RateInter,
2640 fref/MAX2(1e-37,fout),
2641 stage/MAX2(1e-37,fout),
2642 /* simple H+ */
2643 safe_div( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*dense.xIonDense[ipHYDROGEN][0], dense.eden*dense.xIonDense[ipHYDROGEN][1] ),
2644 secondaries.csupra[ipHYDROGEN][0]);
2645
2646 GammaPrt(iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipOpac,
2647 save.ipPnunit[ipPun],iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*
2648 0.05);
2649 }
2650 }
2651
2652 else if( strcmp(save.chSave[ipPun],"HYDp") == 0 )
2653 {
2654 if( strcmp(chTime,"LAST") != 0 )
2655 {
2656 /* save hydrogen populations
2657 * first give total atom and ion density [cm-3]*/
2658 fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e",
2659 radius.depth_mid_zone,
2660 dense.xIonDense[ipHYDROGEN][0],
2661 dense.xIonDense[ipHYDROGEN][1] );
2662
2663 /* next give state-specific densities [cm-3] */
2664 for( j=ipH1s; j < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local-1; j++ )
2665 {
2666 fprintf( save.ipPnunit[ipPun], "\t%.2e",
2667 iso_sp[ipH_LIKE][ipHYDROGEN].st[j].Pop() );
2668 }
2669 fprintf( save.ipPnunit[ipPun], "\n" );
2670 }
2671 }
2672
2673 else if( strcmp(save.chSave[ipPun],"HYDl") == 0 )
2674 {
2675 if( strcmp(chTime,"LAST") == 0 )
2676 {
2677 /* save hydrogen line
2678 * gives intensities and optical depths */
2679 for( ipHi=1; ipHi<iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local -
2680 iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local; ++ipHi )
2681 {
2682 for( ipLo=0; ipLo<ipHi; ++ipLo )
2683 {
2684 if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).ipCont() < 0 )
2685 continue;
2686 fprintf(save.ipPnunit[ipPun], "%li\t%li\t%li\t%li\t%.4e\t%.2e\n",
2687 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].n(),
2688 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].l(),
2689 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].n(),
2690 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].l(),
2691 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).EnergyRyd(),
2692 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).Emis().TauIn() );
2693 }
2694 }
2695 }
2696 }
2697
2698 /* save hydrogen Lya - some details about Lya */
2699 else if( strcmp(save.chSave[ipPun],"HYDL") == 0 )
2700 {
2701 if( strcmp(chTime,"LAST") != 0 )
2702 {
2703 /* the population ratio for Lya */
2704 double popul = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()/SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop());
2705 /* the excitation temperature of Lya */
2706 texc = TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) );
2707 fprintf( save.ipPnunit[ipPun],
2708 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2709 radius.depth_mid_zone,
2710 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2711 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot(),
2712 popul,
2713 texc,
2714 phycon.te,
2715 texc/phycon.te ,
2716 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc(),
2717 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest(),
2718 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().pump(),
2719 opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1],
2720 opac.albedo[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1] );
2721 }
2722 }
2723
2724 else if( strcmp(save.chSave[ipPun],"HYDr") == 0 )
2725 {
2726 /* save hydrogen recc - recombination cooling for AGN3 */
2727 TempChange(2500.f, false);
2728 while( phycon.te <= 20000. )
2729 {
2730 double r1;
2731 double ThinCoolingCaseB;
2732
2733 r1 = HydroRecCool(1,0);
2734 ThinCoolingCaseB = pow(10.,((-25.859117 +
2735 0.16229407*phycon.telogn[0] +
2736 0.34912863*phycon.telogn[1] -
2737 0.10615964*phycon.telogn[2])/(1. +
2738 0.050866793*phycon.telogn[0] -
2739 0.014118924*phycon.telogn[1] +
2740 0.0044980897*phycon.telogn[2] +
2741 6.0969594e-5*phycon.telogn[3])))/phycon.te;
2742
2743 fprintf( save.ipPnunit[ipPun], " %10.2e\t",
2744 phycon.te);
2745 fprintf( save.ipPnunit[ipPun], " %10.2e\t",
2746 (r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) );
2747
2748 fprintf( save.ipPnunit[ipPun], " %10.2e\t",
2749 r1/(BOLTZMANN*phycon.te));
2750
2751 fprintf( save.ipPnunit[ipPun], " %10.2e\n",
2752 ThinCoolingCaseB/(BOLTZMANN*phycon.te));
2753
2754 TempChange(phycon.te *2.f , false);
2755 }
2756 /* must exit since we have disturbed the solution */
2757 fprintf(ioQQQ , "save agn now exits since solution is disturbed.\n");
2758 cdEXIT( EXIT_SUCCESS );
2759 }
2760 else
2761 TotalInsanity();
2762 }
2763
2764 else if( strcmp(save.chSave[ipPun],"IONI") == 0 )
2765 {
2766 if( strcmp(chTime,"LAST") == 0 )
2767 {
2768 /* save mean ionization distribution */
2769 PrtMeanIon( 'i', false , save.ipPnunit[ipPun] );
2770 }
2771 }
2772
2773 /* save ionization rates */
2774 else if( strcmp(save.chSave[ipPun],"IONR") == 0 )
2775 {
2776 if( strcmp(chTime,"LAST") != 0 )
2777 {
2778 /* this is element number */
2779 nelem = (long)save.punarg[ipPun][0];
2780 fprintf( save.ipPnunit[ipPun],
2781 "%.5e\t%.4e\t%.4e",
2782 radius.depth_mid_zone,
2783 dense.eden ,
2784 dynamics.Rate);
2785 /* >>chng 04 oct 15, from nelem+2 to nelem+1 - array over read -
2786 * caught by PnH */
2787 for( ion=0; ion<nelem+1; ++ion )
2788 {
2789 fprintf( save.ipPnunit[ipPun],
2790 "\t%.4e\t%.4e\t%.4e\t%.4e",
2791 dense.xIonDense[nelem][ion] ,
2792 ionbal.RateIonizTot(nelem,ion) ,
2793 ionbal.RateRecomTot[nelem][ion] ,
2794 dynamics.Source[nelem][ion] );
2795 }
2796 fprintf( save.ipPnunit[ipPun], "\n");
2797 }
2798 }
2799
2800 else if( strcmp(save.chSave[ipPun]," IP ") == 0 )
2801 {
2802 if( strcmp(chTime,"LAST") == 0 )
2803 {
2804 /* save valence shell ip's */
2805 for( nelem=0; nelem < LIMELM; nelem++ )
2806 {
2807 int ion_big;
2808 double energy;
2809
2810 /* this is the largest number of ion stages per line */
2811 const int NELEM_LINE = 10;
2812 /* this loop in case all ions do not fit across page */
2813 for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE )
2814 {
2815 int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem);
2816
2817 /* new line then element name */
2818 fprintf( save.ipPnunit[ipPun],
2819 "\n%2.2s", elementnames.chElementSym[nelem]);
2820
2821 /* print ion stages across line */
2822 for( ion=ion_big; ion <= ion_limit; ++ion )
2823 {
2824 fprintf( save.ipPnunit[ipPun], "\t%4ld", ion+1 );
2825 }
2826 fprintf( save.ipPnunit[ipPun], "\n" );
2827
2828 /* this loop is over all shells */
2829 ASSERT( ion_limit < LIMELM );
2830 /* upper limit is number of shells in atom */
2831 for( ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ )
2832 {
2833
2834 /* print shell label */
2835 fprintf( save.ipPnunit[ipPun], "%2.2s", Heavy.chShell[ips]);
2836
2837 /* loop over possible ions */
2838 for( ion=ion_big; ion<=ion_limit; ++ion )
2839 {
2840
2841 /* does this subshell exist for this ion? break if it does not*/
2842 /*if( Heavy.nsShells[nelem][ion]<Heavy.nsShells[nelem][0] )*/
2843 if( ips >= Heavy.nsShells[nelem][ion] )
2844 break;
2845
2846 /* array elements are shell, numb of electrons, element, 0 */
2847 energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0);
2848
2849 /* now print threshold with correct format */
2850 if( energy < 10. )
2851 {
2852 fprintf( save.ipPnunit[ipPun], "\t%6.3f", energy );
2853 }
2854 else if( energy < 100. )
2855 {
2856 fprintf( save.ipPnunit[ipPun], "\t%6.2f", energy );
2857 }
2858 else if( energy < 1000. )
2859 {
2860 fprintf( save.ipPnunit[ipPun], "\t%6.1f", energy );
2861 }
2862 else
2863 {
2864 fprintf( save.ipPnunit[ipPun], "\t%6ld", (long)(energy) );
2865 }
2866 }
2867
2868 /* put cs at end of long line */
2869 fprintf( save.ipPnunit[ipPun], "\n" );
2870 }
2871 }
2872 }
2873 }
2874 }
2875
2876 else if( strcmp(save.chSave[ipPun],"LINC") == 0 )
2877 {
2878 /* save line cumulative */
2879 if( strcmp(chTime,"LAST") != 0 )
2880 {
2881 save_line(save.ipPnunit[ipPun],"PUNC",
2882 save.lgEmergent[ipPun]);
2883 }
2884 }
2885
2886 else if( strcmp(save.chSave[ipPun],"LIND") == 0 )
2887 {
2888 /* save line data, then stop */
2889 SaveLineData(save.ipPnunit[ipPun]);
2890 }
2891
2892 else if( strcmp(save.chSave[ipPun],"LINL") == 0 )
2893 {
2894 /* save line labels */
2895 bool lgPrintAll=false;
2896 /* LONG keyword on save line labels command sets this to 1 */
2897 if( save.punarg[ipPun][0]>0. )
2898 lgPrintAll = true;
2899 prt_LineLabels(save.ipPnunit[ipPun] , lgPrintAll );
2900 }
2901
2902 else if( strcmp(save.chSave[ipPun],"LINO") == 0 )
2903 {
2904 if( strcmp(chTime,"LAST") == 0 )
2905 {
2906 /* save line optical depths, routine is below, file static */
2907 SaveLineStuff(save.ipPnunit[ipPun],"optical" , save.punarg[ipPun][0]);
2908 }
2909 }
2910
2911 else if( strcmp(save.chSave[ipPun],"LINP") == 0 )
2912 {
2913 if( strcmp(chTime,"LAST") != 0 )
2914 {
2915 static bool lgFirst=true;
2916 /* save line populations, need to do this twice if very first
2917 * call since first call to SaveLineStuff generates atomic parameters
2918 * rather than level pops, routine is below, file static */
2919 SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]);
2920 if( lgFirst )
2921 {
2922 lgFirst = false;
2923 SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]);
2924 }
2925 }
2926 }
2927
2928 else if( strcmp(save.chSave[ipPun],"LINS") == 0 )
2929 {
2930 /* save line emissivity */
2931 if( strcmp(chTime,"LAST") != 0 )
2932 {
2933 save_line(save.ipPnunit[ipPun],"PUNS",
2934 save.lgEmergent[ipPun]);
2935 }
2936 }
2937
2938 else if( strcmp(save.chSave[ipPun],"LINR") == 0 )
2939 {
2940 /* save line RT */
2941 if( strcmp(chTime,"LAST") != 0 )
2942 Save_Line_RT( save.ipPnunit[ipPun]);
2943 }
2944
2945 else if( strcmp(save.chSave[ipPun],"LINA") == 0 )
2946 {
2947 /* save line array */
2948 if( strcmp(chTime,"LAST") == 0 )
2949 {
2950 /* save out all lines with energies */
2951 for( j=0; j < LineSave.nsum; j++ )
2952 {
2953 if( LineSv[j].wavelength > 0. &&
2954 LineSv[j].SumLine[0] > 0. )
2955 {
2956 /* line energy, in units set with units option */
2957 fprintf( save.ipPnunit[ipPun], "%12.5e",
2958 AnuUnit((realnum)RYDLAM/LineSv[j].wavelength) );
2959 /* line label */
2960 fprintf( save.ipPnunit[ipPun], "\t%4.4s ",
2961 LineSv[j].chALab );
2962 /* wavelength */
2963 prt_wl( save.ipPnunit[ipPun], LineSv[j].wavelength );
2964 /* intrinsic intensity */
2965 fprintf( save.ipPnunit[ipPun], "\t%8.3f",
2966 log10(SDIV(LineSv[j].SumLine[0]) ) + radius.Conv2PrtInten );
2967 /* emergent line intensity, r recombination */
2968 fprintf( save.ipPnunit[ipPun], "\t%8.3f",
2969 log10(SDIV(LineSv[j].SumLine[1]) ) + radius.Conv2PrtInten );
2970 /* type of line, i for info, etc */
2971 fprintf( save.ipPnunit[ipPun], " \t%c\n",
2972 LineSv[j].chSumTyp);
2973 }
2974 }
2975 }
2976 }
2977
2978 else if( strcmp(save.chSave[ipPun],"LINI") == 0 )
2979 {
2980 if( strcmp(chTime,"LAST") == 0 &&
2981 (nzone/save.LinEvery)*save.LinEvery != nzone )
2982 {
2983 /* this is the last zone
2984 * save line intensities - but do not do last zone twice */
2985 SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0] );
2986 }
2987 else if( strcmp(chTime,"LAST") != 0 )
2988 {
2989 /* following so we only save first zone if LinEvery reset */
2990 if( (save.lgLinEvery && nzone == 1) ||
2991 (nzone/save.LinEvery)*save.LinEvery == nzone )
2992 {
2993 /* this is middle of calculation
2994 * save line intensities */
2995 SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0]);
2996 }
2997 }
2998 }
2999
3000 else if( strcmp( save.chSave[ipPun],"LEIL") == 0)
3001 {
3002 /* some line intensities for the Leiden PDR,
3003 * but only do this when calculation is complete */
3004 if( strcmp(chTime,"LAST") == 0 )
3005 {
3006 double absval , rel;
3007 long int n;
3008 /* the lines we will find,
3009 * for a sample list of PDR lines look at LineList_PDR_H2.dat
3010 * in the cloudy data dir */
3011 /* the number of H2 lines */
3012 const int NLINE_H2 = 31;
3013 /* the number of lines which are not H2 */
3014 const int NLINE_NOTH_H2 = 5;
3015 /* the labels and wavelengths for the lines that are not H2 */
3016 char chLabel[NLINE_NOTH_H2][5]=
3017 {"C 2", "O 1", "O 1","C 1", "C 1" };
3018 double Wl[NLINE_NOTH_H2]=
3019 {157.6 , 63.17 , 145.5 ,609.2 , 369.7 , };
3020 /* these are wavelengths in microns, conv to Angstroms before call */
3021 /* >>chng 05 sep 06, many of following wavelengths updated to agree
3022 * with output - apparently not updated when energies changed */
3023 double Wl_H2[NLINE_H2]=
3024 {2.121 ,
3025 28.213, 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 ,
3026 4.408 , 4.180 , 3.996 , 3.845 , 3.723 , 3.625 , 3.546 , 3.485 , 3.437 , 3.403 ,
3027 3.380 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603};
3028 /* print a header for the lines */
3029 for( n=0; n<NLINE_NOTH_H2; ++n )
3030 {
3031 fprintf(save.ipPnunit[ipPun], "%s\t%.2f",chLabel[n] , Wl[n]);
3032 /* get the line, non positive return says didn't find it */
3033 /* arguments are 4-char label, wavelength, return log total intensity, linear rel inten */
3034 if( cdLine( chLabel[n] , (realnum)(Wl[n]*1e4) , &absval , &rel ) <= 0 )
3035 {
3036 fprintf(save.ipPnunit[ipPun], " did not find\n");
3037 }
3038 else
3039 {
3040 fprintf(save.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
3041 }
3042 }
3043 fprintf(save.ipPnunit[ipPun], "\n\n\n");
3044
3045 /* only print the H2 lines if the big molecule is turned on */
3046 if( h2.lgEnabled )
3047 {
3048 fprintf(save.ipPnunit[ipPun],
3049 "Here are some of the H2 Intensities, The first one is the\n"
3050 "1-0 S(0) line and the following ones are the 0-0 S(X)\n"
3051 "lines where X goes from 0 to 29\n\n");
3052 for( n=0; n<NLINE_H2; ++n )
3053 {
3054 fprintf(save.ipPnunit[ipPun], "%s\t%.3f","H2 " , Wl_H2[n]);
3055 /* get the line, non positive return says didn't find it */
3056 if( cdLine( "H2 " , (realnum)(Wl_H2[n]*1e4) , &absval , &rel ) <= 0 )
3057 {
3058 fprintf(save.ipPnunit[ipPun], " did not find\n");
3059 }
3060 else
3061 {
3062 fprintf(save.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
3063 }
3064 }
3065 }
3066 }
3067 }
3068
3069 else if( strcmp( save.chSave[ipPun],"LEIS") == 0)
3070 {
3071 if( strcmp(chTime,"LAST") != 0 )
3072 {
3073 /* get some column densities we shall need */
3074 double col_ci , col_oi , col_cii, col_heii;
3075 if( cdColm("carb" , 1 , &col_ci ) )
3076 TotalInsanity();
3077 if( cdColm("carb" , 2 , &col_cii ) )
3078 TotalInsanity();
3079 if( cdColm("oxyg" , 1 , &col_oi ) )
3080 TotalInsanity();
3081 if( cdColm("heli" , 2 , &col_heii ) )
3082 TotalInsanity();
3083 /* save Leiden structure - some numbers for the Leiden PDR model comparisons */
3084 fprintf( save.ipPnunit[ipPun],
3085 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
3086 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
3087 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
3088 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
3089 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
3090 /* depth of this point */
3091 radius.depth_mid_zone,
3092 /* A_V for an extended source */
3093 0.00,
3094 /* A_V for a point source */
3095 rfield.extin_mag_V_point,
3096 /* temperature */
3097 phycon.te ,
3098 dense.xIonDense[ipHYDROGEN][0],
3099 hmi.H2_total,
3100 dense.xIonDense[ipCARBON][0],
3101 dense.xIonDense[ipCARBON][1],
3102 dense.xIonDense[ipOXYGEN][0],
3103 findspecieslocal("CO")->den,
3104 findspecieslocal("O2")->den,
3105 findspecieslocal("CH")->den,
3106 findspecieslocal("OH")->den,
3107 dense.eden,
3108 dense.xIonDense[ipHELIUM][1],
3109 dense.xIonDense[ipHYDROGEN][1],
3110 findspecieslocal("H3+")->den,
3111 colden.colden[ipCOL_H0],
3112 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
3113 col_ci,
3114 col_cii,
3115 col_oi,
3116 findspecieslocal("CO")->column,
3117 findspecieslocal("O2")->column,
3118 findspecieslocal("CH")->column,
3119 findspecieslocal("OH")->column,
3120 colden.colden[ipCOL_elec],
3121 col_heii,
3122 colden.colden[ipCOL_Hp],
3123 colden.colden[ipCOL_H3p],
3124 hmi.H2_Solomon_dissoc_rate_used_H2g ,
3125 gv.rate_h2_form_grains_used_total,
3126 hmi.H2_photodissoc_used_H2g,
3127 hmi.UV_Cont_rel2_Draine_DB96_depth,
3128 /* CO and C dissociation rate */
3129 mole.findrk("PHOTON,CO=>C,O"),
3130 /* total CI ionization rate */
3131 ionbal.PhotoRate_Shell[ipCARBON][0][2][0],
3132 /* total heating, erg cm-3 s-1 */
3133 thermal.htot,
3134 /* total cooling, erg cm-3 s-1 */
3135 thermal.ctot,
3136 /* GrnP grain photo heating */
3137 thermal.heating[0][13],
3138 /* grain collisional cooling */
3139 MAX2(0.,gv.GasCoolColl),
3140 /* grain collisional heating */
3141 -1.*MIN2(0.,gv.GasCoolColl),
3142 /* COds - CO dissociation heating */
3143 thermal.heating[0][9],
3144 /* H2dH-Heating due to H2 dissociation */
3145 hmi.HeatH2Dish_used,
3146 /* H2vH-Heating due to collisions with H2 */
3147 hmi.HeatH2Dexc_used ,
3148 /* ChaT - charge transfer heating */
3149 thermal.heating[0][24] ,
3150 /* cosmic ray heating */
3151 thermal.heating[1][6] ,
3152 /* heating due to atoms of various heavy elements */
3153 thermal.heating[ipMAGNESIUM][0],
3154 thermal.heating[ipSULPHUR][0],
3155 thermal.heating[ipSILICON][0],
3156 thermal.heating[ipIRON][0],
3157 thermal.heating[ipSODIUM][0],
3158 thermal.heating[ipALUMINIUM][0],
3159 thermal.heating[ipCARBON][0],
3160 TauLines[ipT610].Coll().cool(),
3161 TauLines[ipT370].Coll().cool(),
3162 TauLines[ipT157].Coll().cool(),
3163 TauLines[ipT63].Coll().cool(),
3164 TauLines[ipT146].Coll().cool() );
3165 }
3166 }
3167
3168# ifdef USE_NLTE7
3169 else if( strcmp( save.chSave[ipPun],"NLTE") == 0)
3170 {
3171 if( strcmp(chTime,"LAST") == 0 )
3172 {
3173 //Generate output for NLTE7 conference
3174 runNLTE(ipPun);
3175 }
3176 }
3177# endif
3178
3179
3180 /* Print out the FeII energy level file into Stout format*/
3181 else if( strcmp( save.chSave[ipPun],"LY1") == 0)
3182 {
3183 static bool runonce = true;
3184 if (runonce )
3185 {
3186 for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
3187 {
3188 ASSERT(FeII.FeIINRGs[ipHi] >= 0.);
3189 ASSERT(FeII.FeIISTWT[ipHi] >= 0.);
3190 fprintf(save.ipPnunit[ipPun], "%li\t%10.3f\t%3.1f\n",ipHi+1,FeII.FeIINRGs[ipHi],FeII.FeIISTWT[ipHi]);
3191 }
3192 runonce = false;
3193 fprintf(save.ipPnunit[ipPun], "-1\n");
3194 }
3195 }
3196
3197 /* Print out the FeII transition probablility file into Stout format*/
3198 else if( strcmp( save.chSave[ipPun],"LY2") == 0)
3199 {
3200 static bool runonce = true;
3201 if (runonce )
3202 {
3203 for( long ipLo = 0; ipLo < FeII.nFeIILevel_malloc; ipLo++)
3204 {
3205 for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
3206 {
3207 if( FeII.FeIIAul[ipHi][ipLo] >= 0. && FeII.FeIIAul[ipHi][ipLo] != 1e-5f && FeII.FeIIAul[ipHi][ipLo] != 1e-6f )
3208 {
3209 fprintf(save.ipPnunit[ipPun], "A\t%li\t%li\t%8.2e\n",
3210 ipLo+1,ipHi+1,FeII.FeIIAul[ipHi][ipLo]);
3211 }
3212 }
3213 }
3214 runonce = false;
3215 fprintf(save.ipPnunit[ipPun], "-1\n");
3216 }
3217 }
3218
3219 /* Print out the FeII Collision data file into Stout format*/
3220 else if( strcmp( save.chSave[ipPun],"LY3") == 0)
3221 {
3222 static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
3223 static bool runonce = true;
3224 if (runonce )
3225 {
3226 fprintf(save.ipPnunit[ipPun], "TEMP\t8");
3227 for( int k = 0; k < 8; k++)
3228 {
3229 fprintf(save.ipPnunit[ipPun], "\t%8.2e",tt[k]);
3230 }
3231 fprintf(save.ipPnunit[ipPun], "\n");
3232 for( long ipHi = 1; ipHi < FeII.nFeIILevel_malloc; ipHi++)
3233 {
3234 for( long ipLo=0; ipLo < ipHi; ipLo++ )
3235 {
3236 bool skipLine = false;
3237 for( int k = 0; k < 8; k++)
3238 {
3239 if( FeII.FeIIColl[ipHi][ipLo][k] == -2. || FeII.FeIIColl[ipHi][ipLo][k] == 0.)
3240 {
3241 skipLine = true;
3242 break;
3243 }
3244 }
3245 if( !skipLine )
3246 {
3247 fprintf(save.ipPnunit[ipPun], "CSELECTRON\t%li\t%li",ipLo+1,ipHi+1);
3248 for( int k = 0; k < 8; k++)
3249 {
3250 fprintf(save.ipPnunit[ipPun], "\t%8.2e",FeII.FeIIColl[ipHi][ipLo][k]);
3251 }
3252 fprintf(save.ipPnunit[ipPun], "\n");
3253 }
3254 }
3255 }
3256 runonce = false;
3257 fprintf(save.ipPnunit[ipPun], "-1\n");
3258 }
3259 }
3260
3261 else if( strcmp( save.chSave[ipPun],"LLST") == 0)
3262 {
3263 /* save linelist command - do on last iteration */
3264 if( strcmp(chTime,"LAST") == 0 )
3265 {
3266 fprintf( save.ipPnunit[ipPun], "iteration %li" , iteration );
3267 if( save.punarg[ipPun][1] )// column print
3268 fprintf( save.ipPnunit[ipPun], "\n" );
3269
3270 /* -1 is flag saying that this save command was not set */
3271 if( save.nLineList[ipPun] < 0 )
3272 TotalInsanity();
3273
3274 int LineType = 0;
3275 if( save.lgEmergent[ipPun] )
3276 LineType = 1;
3277 if( save.lgCumulative[ipPun] )
3278 LineType += 2;
3279
3280 /* loop over all lines in the file we read */
3281 for( j=0; j<save.nLineList[ipPun]; ++j )
3282 {
3283 double relative , absolute, PrtQuantity;
3284 if( (cdLine( save.chLineListLabel[ipPun][j] ,
3285 save.wlLineList[ipPun][j] ,
3286 &relative , &absolute , LineType ) ) <=0 )
3287 {
3288 if( !h2.lgEnabled && strncmp( save.chLineListLabel[ipPun][j] , "H2 " , 4 )==0 )
3289 {
3290 static bool lgMustPrintFirstTime = true;
3291 if( lgMustPrintFirstTime )
3292 {
3293 /* it's an H2 line and H2 is not being done - ignore it */
3294 fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
3295 "included, so I will ignore it. Log intensity set to -30.\n" );
3296 fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
3297 lgMustPrintFirstTime = false;
3298 }
3299 relative = -30.f;
3300 absolute = -30.f;
3301 }
3302 else if( lgAbort )
3303 {
3304 /* we are in abort mode */
3305 relative = -30.f;
3306 absolute = -30.f;
3307 }
3308 else
3309 {
3310 fprintf(ioQQQ,"DISASTER - did not find a line in the Line List table\n");
3311 cdEXIT(EXIT_FAILURE);
3312 }
3313 }
3314
3315 /* options to do either relative or absolute intensity
3316 * default is relative, is absolute keyword on line then
3317 * punarg set to 1 */
3318 /* straight line intensities */
3319 if( save.punarg[ipPun][0] > 0 )
3320 PrtQuantity = pow(10. , absolute);
3321 else
3322 PrtQuantity = relative;
3323
3324 // column mode, print label
3325 if( save.punarg[ipPun][1] )
3326 {
3327 /* if taking ratio then put div sign between pairs */
3328 if( save.lgLineListRatio[ipPun] && is_odd(j) )
3329 fprintf( save.ipPnunit[ipPun] , "/" );
3330
3331 fprintf( save.ipPnunit[ipPun], "%s ", save.chLineListLabel[ipPun][j] );
3332 char chTemp[MAX_HEADER_SIZE];
3333 sprt_wl( chTemp, save.wlLineList[ipPun][j] );
3334 fprintf( save.ipPnunit[ipPun], "%s ", chTemp );
3335 }
3336
3337 /* if taking ratio print every other line as ratio
3338 * with previous line */
3339 if( save.lgLineListRatio[ipPun] )
3340 {
3341 /* do line pair ratios */
3342 static double SaveQuantity = 0;
3343 if( is_odd(j) )
3344 fprintf( save.ipPnunit[ipPun], "\t%.4e" ,
3345 SaveQuantity / SDIV( PrtQuantity ) );
3346 else
3347 SaveQuantity = PrtQuantity;
3348 }
3349 else
3350 {
3351 fprintf( save.ipPnunit[ipPun], "\t%.4e" , PrtQuantity );
3352 }
3353 // column printout, but check if first of pair
3354 if( save.punarg[ipPun][1] )
3355 {
3356 if( !save.lgLineListRatio[ipPun] ||
3357 is_odd(j) )
3358 fprintf( save.ipPnunit[ipPun], "\n" );
3359 }
3360 }
3361 fprintf( save.ipPnunit[ipPun], "\n" );
3362 }
3363 }
3364
3365 else if( strcmp( save.chSave[ipPun],"CHRT") == 0)
3366 {
3367 /* save chemistry rates command */
3368 if( strcmp(chTime,"LAST") != 0 )
3369 {
3370 bool lgHeader, lgData;
3371 if( save.lgPunHeader[ipPun] )
3372 {
3373 lgHeader = true;
3374 lgData = false;
3375 mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone);
3376 }
3377 save.lgPunHeader[ipPun] = false;
3378 lgHeader = false;
3379 lgData = true;
3380 mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone);
3381 }
3382 }
3383
3384 else if( strcmp(save.chSave[ipPun],"MAP ") == 0 )
3385 {
3386 /* do the map now if we are at the zone, or if this
3387 * is the LAST call to this routine and map not done yet */
3388 if( !hcmap.lgMapDone &&
3389 (nzone == hcmap.MapZone || strcmp(chTime,"LAST") == 0 ) )
3390 {
3391 lgTlkSav = called.lgTalk;
3392 called.lgTalk = cpu.i().lgMPI_talk();
3393 hcmap.lgMapBeingDone = true;
3394 map_do(save.ipPnunit[ipPun] , " map");
3395 called.lgTalk = lgTlkSav;
3396 }
3397 }
3398
3399 else if( strcmp(save.chSave[ipPun],"MOLE") == 0 )
3400 {
3401 if( save.lgPunHeader[ipPun] )
3402 {
3403 fprintf( save.ipPnunit[ipPun],
3404 "#molecular species will follow:\n");
3405 fprintf( save.ipPnunit[ipPun],
3406 "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate");
3407
3408 for(i=0; i<mole_global.num_calc; ++i )
3409 {
3410 fprintf( save.ipPnunit[ipPun], "\t%s", mole_global.list[i]->label.c_str() );
3411 }
3412 fprintf ( save.ipPnunit[ipPun], "\n");
3413 save.lgPunHeader[ipPun] = false;
3414 }
3415 if( strcmp(chTime,"LAST") != 0 )
3416 {
3417 /* molecules, especially for PDR, first give radius */
3418 fprintf( save.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone );
3419
3420 /* visual extinction of point source (star)*/
3421 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
3422
3423 /* visual extinction of an extended source (like a PDR)*/
3424 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
3425
3426 /* carbon monoxide photodissociation rate */
3427 fprintf( save.ipPnunit[ipPun], "%.5e\t" , mole.findrk("PHOTON,CO=>C,O") );
3428
3429 /* carbon recombination rate */
3430 fprintf( save.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] );
3431
3432 /* now do all the molecules */
3433 for(j=0; j<mole_global.num_calc; ++j )
3434 {
3435 fprintf(save.ipPnunit[ipPun],"\t%.2e",mole.species[j].den );
3436 }
3437
3438 fprintf(save.ipPnunit[ipPun],"\n");
3439 }
3440 }
3441
3442 else if( strcmp(save.chSave[ipPun],"OPAC") == 0 )
3443 {
3444 /* save opacity- routine will parse which type of opacity save to do */
3445 if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
3446 save_opacity(save.ipPnunit[ipPun],ipPun);
3447 }
3448
3449 /* save coarse optical depths command */
3450 else if( strcmp(save.chSave[ipPun],"OPTc") == 0 )
3451 {
3452 if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
3453 {
3454 for( j=0; j < rfield.nflux; j++ )
3455 {
3456 fprintf( save.ipPnunit[ipPun],
3457 "%13.5e\t%.3e\t%12.4e\t%.3e\n",
3458 AnuUnit(rfield.AnuOrg[j]),
3459 opac.TauAbsFace[j]+opac.TauScatFace[j],
3460 opac.TauAbsFace[j],
3461 opac.TauScatFace[j] );
3462 }
3463 }
3464 }
3465
3466 /* save fine optical depths command */
3467 else if( strcmp(save.chSave[ipPun],"OPTf") == 0 )
3468 {
3469 if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
3470 {
3471 long nu_hi , nskip;
3472 if( save.punarg[ipPun][0] > 0. )
3473 /* possible lower bounds to energy range - will be zero if not set */
3474 j = ipFineCont( save.punarg[ipPun][0] );
3475 else
3476 j = 0;
3477
3478 /* upper limit */
3479 if( save.punarg[ipPun][1]> 0. )
3480 nu_hi = ipFineCont( save.punarg[ipPun][1]);
3481 else
3482 nu_hi = rfield.nfine;
3483
3484 /* we will bring nskip cells together into one printed
3485 * number to make output smaller - default is 10 */
3486 nskip = (long)save.punarg[ipPun][2];
3487 do
3488 {
3489 realnum sum1 = rfield.fine_opt_depth[j];
3490 realnum sum2 = rfield.fine_opac_zone[j];
3491 /* want to report the central wavelength of the cell */
3492 realnum xnu = rfield.fine_anu[j];
3493 for( jj=1; jj<nskip; ++jj )
3494 {
3495 sum1 += rfield.fine_opt_depth[j+jj];
3496 sum2 += rfield.fine_opac_zone[j+jj];
3497 xnu += rfield.fine_anu[j+jj];
3498 }
3499 if( sum2 > 0. )
3500 fprintf( save.ipPnunit[ipPun],
3501 "%12.6e\t%.3e\t%.3e\n",
3502 AnuUnit(xnu/nskip),
3503 sum1/nskip ,
3504 sum2/nskip);
3505 j += nskip;
3506 }while( j < nu_hi );
3507 }
3508 }
3509
3510 else if( strcmp(save.chSave[ipPun]," OTS") == 0 )
3511 {
3512 ConMax = 0.;
3513 xLinMax = 0.;
3514 opConSum = 0.;
3515 opLinSum = 0.;
3516 ipLinMax = 1;
3517 ipConMax = 1;
3518
3519 for( j=0; j < rfield.nflux; j++ )
3520 {
3521 opConSum += rfield.otscon[j]*opac.opacity_abs[j];
3522 opLinSum += rfield.otslin[j]*opac.opacity_abs[j];
3523 if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax )
3524 {
3525 xLinMax = rfield.otslin[j]*opac.opacity_abs[j];
3526 ipLinMax = j+1;
3527 }
3528 if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax )
3529 {
3530 ConMax = rfield.otscon[j]*opac.opacity_abs[j];
3531 ipConMax = j+1;
3532 }
3533 }
3534 fprintf( save.ipPnunit[ipPun],
3535 "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n",
3536 opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1]
3537 , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1]
3538 , rfield.anu[ipConMax-1], ConMax );
3539 }
3540
3541 else if( strcmp(save.chSave[ipPun],"OVER") == 0 )
3542 {
3543 /* save overview
3544 * this is the floor for the smallest ionization fractions printed */
3545 double toosmall = SMALLFLOAT ,
3546 hold;
3547
3548 /* overview of model results,
3549 * depth, te, hden, eden, ion fractions H, He, c, O */
3550 if( strcmp(chTime,"LAST") != 0 )
3551 {
3552
3553 /* print the depth */
3554 fprintf( save.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone );
3555
3556 /* temperature, heating */
3557 if(dynamics.Cool() > dynamics.Heat())
3558 {
3559 fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f",
3560 log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat())) );
3561 }
3562 else
3563 {
3564 double diff = fabs(thermal.htot-dynamics.Cool());
3565 fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f",
3566 log10(phycon.te), log10(SDIV(diff)) );
3567 }
3568
3569 /* hydrogen and electron densities */
3570 fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f",
3571 log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) );
3572
3573 /* molecular fraction of hydrogen */
3574 fprintf( save.ipPnunit[ipPun], "\t%.4f",
3575 /*log10(MAX2(toosmall,2.*findspecies("H2")->den/dense.gas_phase[ipHYDROGEN])) );*/
3576 log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) );
3577
3578 /* ionization fractions of hydrogen */
3579 fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f",
3580 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])),
3581 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );
3582
3583 /* ionization fractions of helium */
3584 for( j=1; j <= 3; j++ )
3585 {
3586 double arg1 = SDIV(dense.gas_phase[ipHELIUM]);
3587 arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 );
3588 fprintf( save.ipPnunit[ipPun], "\t%.4f",
3589 log10(arg1) );
3590 }
3591
3592 /* carbon monoxide molecular fraction of CO */
3593 hold = SDIV(dense.gas_phase[ipCARBON]);
3594 hold = findspecieslocal("CO")->den/hold;
3595 hold = MAX2(toosmall, hold );
3596 fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) );
3597
3598 /* ionization fractions of carbon */
3599 for( j=1; j <= 4; j++ )
3600 {
3601 hold = SDIV(dense.gas_phase[ipCARBON]);
3602 hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold);
3603 fprintf( save.ipPnunit[ipPun], "\t%.4f",
3604 log10(hold) );
3605 }
3606
3607 /* ionization fractions of oxygen */
3608 for( j=1; j <= 6; j++ )
3609 {
3610 hold = SDIV(dense.gas_phase[ipOXYGEN]);
3611 hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold);
3612 fprintf( save.ipPnunit[ipPun], "\t%.4f",
3613 log10(hold) );
3614 }
3615
3616 // molecular fraction of H2O
3617 hold = SDIV(dense.gas_phase[ipOXYGEN]);
3618 hold = findspecieslocal("H2O")->den/hold;
3619 hold = MAX2(toosmall, hold );
3620 fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) );
3621
3622 /* visual extinction of point source (star)*/
3623 fprintf( save.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point);
3624
3625 /* visual extinction of an extended source (like a PDR)*/
3626 fprintf( save.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended);
3627 }
3628 }
3629
3630 else if( strcmp(save.chSave[ipPun]," PDR") == 0 )
3631 {
3632 /* this is the save PDR command */
3633 if( strcmp(chTime,"LAST") != 0 )
3634 {
3635 /* convert optical depth at wavelength of V filter
3636 * into magnitudes of extinction */
3637 /* >>chyng 03 feb 25, report extinction to illuminated face,
3638 * rather than total extinction which included far side when
3639 * sphere was set */
3640 /*av = opac.TauTotalGeo[0][rfield.ipV_filter-1]*1.08574;*/
3641
3642 fprintf( save.ipPnunit[ipPun],
3643 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t",
3644 radius.depth_mid_zone,
3645 /* total hydrogen column density, all forms */
3646 colden.colden[ipCOL_HTOT],
3647 phycon.te,
3648 /* fraction of H that is atomic */
3649 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
3650 /* ratio of n(H2) to total H, == 0.5 when fully molecular */
3651 2.*findspecieslocal("H2")->den/dense.gas_phase[ipHYDROGEN],
3652 2.*findspecieslocal("H2*")->den/dense.gas_phase[ipHYDROGEN],
3653 /* atomic to total carbon */
3654 dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]),
3655 findspecieslocal("CO")->den/SDIV(dense.gas_phase[ipCARBON]),
3656 findspecieslocal("H2O")->den/SDIV(dense.gas_phase[ipOXYGEN]),
3657 /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */
3658 hmi.UV_Cont_rel2_Habing_TH85_depth);
3659
3660 /* visual extinction due to dust alone, of point source (star)*/
3661 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
3662
3663 /* visual extinction due to dust alone, of an extended source (like a PDR)*/
3664 fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
3665
3666 /* visual extinction (all sources) of a point source (like a PDR)*/
3667 fprintf( save.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] );
3668 }
3669 }
3670
3671 /* performance characteristics per zone */
3672 else if( strcmp(save.chSave[ipPun],"PERF") == 0 )
3673 {
3674 if( strcmp(chTime,"LAST") != 0 )
3675 {
3676 static double ElapsedTime , ZoneTime;
3677 if( nzone<=1 )
3678 {
3679 ElapsedTime = cdExecTime();
3680 ZoneTime = 0.;
3681 }
3682 else
3683 {
3684 double t = cdExecTime();
3685 ZoneTime = t - ElapsedTime;
3686 ElapsedTime = t;
3687 }
3688
3689 /* zone, time for this zone, elapsed time */
3690 fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\t%li",
3691 nzone, ZoneTime , ElapsedTime, conv.nPres2Ioniz );
3692 // print various loop counters
3693 for( long i=0; i<NTYPES; ++i )
3694 fprintf( save.ipPnunit[ipPun], "\t%li", conv.getCounterZone(i) );
3695 fprintf( save.ipPnunit[ipPun], "\n" );
3696 }
3697 }
3698
3699 else if( strcmp(save.chSave[ipPun],"PHYS") == 0 )
3700 {
3701 if( strcmp(chTime,"LAST") != 0 )
3702 {
3703 /* save physical conditions */
3704 fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
3705 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN],
3706 dense.eden, thermal.htot, wind.AccelTotalOutward, geometry.FillFac );
3707 }
3708 }
3709
3710 else if( strcmp(save.chSave[ipPun],"PRES") == 0 )
3711 {
3712 /* the save pressure command */
3713 if( strcmp(chTime,"LAST") != 0 )
3714 {
3715 fprintf( save.ipPnunit[ipPun],
3716 "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%c\n",
3717 /*A 1 #P depth */
3718 radius.depth_mid_zone,
3719 /*B 2 Perror */
3720 pressure.PresTotlError*100.,
3721 /*C 3 Pcurrent */
3722 pressure.PresTotlCurr,
3723 /*D 4 Pln + pintg
3724 * >>chng 06 apr 19, subtract pinzon the acceleration added in this zone
3725 * since is not total at outer edge of zone, above is at inner edge */
3726 pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon,
3727 /*E 5 pgas (0) */
3728 pressure.PresTotlInit,
3729 /*F 6 Pgas */
3730 pressure.PresGasCurr,
3731 /*G 7 Pram */
3732 pressure.PresRamCurr,
3733 /*H 8 P rad in lines */
3734 pressure.pres_radiation_lines_curr,
3735 /*I 9 Pinteg subtract continuum rad pres which has already been added on */
3736 pressure.PresInteg - pressure.pinzon,
3737 /*J 10 V(wind km/s) wind speed in km/s */
3738 wind.windv/1e5,
3739 /*K cad(km/s) sound speed in km/s */
3740 timesc.sound_speed_adiabatic/1e5,
3741 /* the magnetic pressure */
3742 magnetic.pressure ,
3743 /* the local turbulent velocity in km/s */
3744 DoppVel.TurbVel/1e5 ,
3745 /* turbulent pressure */
3746 pressure.PresTurbCurr*DoppVel.lgTurb_pressure ,
3747 /* gravitational pressure */
3748 pressure.IntegRhoGravity,
3749 // the integral of electron scattering acceleration in
3750 // the absence of any absorptio, minus acceleration in current
3751 // zone, which has been added in - done this way, result is
3752 // zero in first zonen
3753 pressure.PresIntegElecThin-pressure.pinzon_PresIntegElecThin,
3754 // is this converged?
3755 TorF(conv.lgConvPres) );
3756 }
3757 }
3758 else if( strcmp(save.chSave[ipPun],"PREL") == 0 )
3759 {
3760 /* line pressure contributors */
3761 fprintf( save.ipPnunit[ipPun],
3762 "%.5e\t%.3e\t%.3e\t",
3763 /*A 1 #P depth */
3764 radius.depth_mid_zone ,
3765 pressure.PresTotlCurr,
3766 pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr) );
3767 PrtLinePres(save.ipPnunit[ipPun]);
3768
3769 }
3770
3771 else if( save.chSave[ipPun][0]=='R' )
3772 {
3773 /* work around internal limits to Microsoft vs compiler */
3774 if( strcmp(save.chSave[ipPun],"RADI") == 0 )
3775 {
3776 /* save radius information for all zones */
3777 if( strcmp(chTime,"LAST") != 0 )
3778 {
3779 fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
3780 nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
3781 radius.drad );
3782 }
3783 }
3784
3785 else if( strcmp(save.chSave[ipPun],"RADO") == 0 )
3786 {
3787 /* save radius information for only the last zone */
3788 if( strcmp(chTime,"LAST") == 0 )
3789 {
3790 fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
3791 nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
3792 radius.drad );
3793 }
3794 }
3795
3796 else if( strcmp(save.chSave[ipPun],"RESU") == 0 )
3797 {
3798 /* save results of the calculation */
3799 if( strcmp(chTime,"LAST") == 0 )
3800 SaveResults(save.ipPnunit[ipPun]);
3801 }
3802 else
3803 {
3804 /* this can't happen */
3805 TotalInsanity();
3806 }
3807 }
3808
3809 else if( strcmp(save.chSave[ipPun],"SECO") == 0 )
3810 {
3811 /* save secondary ionization */
3812 if( strcmp(chTime,"LAST") != 0 )
3813 fprintf(save.ipPnunit[ipPun],
3814 "%.5e\t%.3e\t%.3e\t%.3e\n",
3815 radius.depth ,
3816 secondaries.csupra[ipHYDROGEN][0],
3817 secondaries.csupra[ipHYDROGEN][0]*2.02,
3818 secondaries.x12tot );
3819 }
3820
3821 else if( strcmp(save.chSave[ipPun],"SOUS") == 0 )
3822 {
3823 /* full spectrum of continuum source function at 1 depth
3824 * command was "save source spectrum" */
3825 if( strcmp(chTime,"LAST") != 0 )
3826 {
3827 limit = MIN2(rfield.ipMaxBolt,rfield.nflux);
3828 for( j=0; j < limit; j++ )
3829 {
3830 fprintf( save.ipPnunit[ipPun],
3831 "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
3832 AnuUnit(rfield.AnuOrg[j]),
3833 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j],
3834 opac.opacity_abs[j],
3835 rfield.ConSourceFcnLocal[nzone][j],
3836 rfield.ConSourceFcnLocal[nzone][j]/plankf(j) ,
3837 safe_div(rfield.ConSourceFcnLocal[nzone][j],rfield.flux[0][j]));
3838 }
3839 }
3840 }
3841
3842 else if( strcmp(save.chSave[ipPun],"SOUD") == 0 )
3843 {
3844 /* parts of continuum source function vs depth
3845 * command was save source function depth */
3846 j = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon + 2;
3847 fprintf( save.ipPnunit[ipPun],
3848 "%.4e\t%.4e\t%.4e\t%.4e\n",
3849 opac.TauAbsFace[j-1],
3850 rfield.ConEmitLocal[nzone][j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]),
3851 rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1],
3852 rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] );
3853 }
3854
3855 /* this is save special option */
3856 else if( strcmp(save.chSave[ipPun],"SPEC") == 0 )
3857 {
3858 SaveSpecial(save.ipPnunit[ipPun],chTime);
3859 }
3860
3861 /* this is save species option */
3862 else if( strcmp(save.chSave[ipPun],"SPCS") == 0 )
3863 {
3864 if( ( strcmp(chTime,"LAST") != 0 && strcmp(save.chSaveArgs[ipPun],"COLU") != 0 ) ||
3865 ( strcmp(chTime,"LAST") == 0 && strcmp(save.chSaveArgs[ipPun],"COLU") == 0 ) )
3866 SaveSpecies(save.ipPnunit[ipPun] , ipPun);
3867 }
3868
3869 else if( strcmp(save.chSave[ipPun],"TEMP") == 0 )
3870 {
3871 static double deriv_old=-1;
3872 double deriv=-1. , deriv_sec;
3873 /* temperature and its derivatives */
3874 fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e",
3875 radius.depth_mid_zone,
3876 phycon.te,
3877 thermal.dCooldT );
3878 /* if second zone then have one deriv */
3879 if( nzone >1 )
3880 {
3881 deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
3882 fprintf( save.ipPnunit[ipPun], "\t%.2e", deriv );
3883 /* if third zone then have second deriv */
3884 if( nzone > 2 )
3885 {
3886 deriv_sec = (deriv-deriv_old)/ radius.drad;
3887 fprintf( save.ipPnunit[ipPun], "\t%.2e",
3888 deriv_sec );
3889 }
3890 deriv_old = deriv;
3891 }
3892 fprintf( save.ipPnunit[ipPun], "\n");
3893 }
3894
3895 /* time dependent model */
3896 else if( strcmp(save.chSave[ipPun],"TIMD") == 0 )
3897 {
3898 if( strcmp(chTime,"LAST") == 0 )
3899 DynaPunchTimeDep( save.ipPnunit[ipPun] , "END" );
3900 }
3901
3902 /* execution time per zone */
3903 else if( strcmp(save.chSave[ipPun],"XTIM") == 0 )
3904 {
3905 static double ElapsedTime , ZoneTime;
3906 if( nzone<=1 )
3907 {
3908 ElapsedTime = cdExecTime();
3909 ZoneTime = 0.;
3910 }
3911 else
3912 {
3913 double t = cdExecTime();
3914 ZoneTime = t - ElapsedTime;
3915 ElapsedTime = t;
3916 }
3917
3918 /* zone, time for this zone, elapsed time */
3919 fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n",
3920 nzone, ZoneTime , ElapsedTime );
3921 }
3922
3923 else if( strcmp(save.chSave[ipPun],"TPRE") == 0 )
3924 {
3925 /* temperature and its predictors, turned on with save tprid */
3926 fprintf( save.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n",
3927 nzone, phycon.TeInit, phycon.TeProp, phycon.te,
3928 (phycon.TeProp- phycon.te)/phycon.te );
3929 }
3930
3931 else if( strcmp(save.chSave[ipPun],"WIND") == 0 )
3932 {
3933 /* wind velocity, radiative acceleration, and ratio total
3934 * to electron scattering acceleration */
3935 /* first test only save last zone */
3936 if( (save.punarg[ipPun][0] == 0 && strcmp(chTime,"LAST") == 0)
3937 ||
3938 /* this test save all zones */
3939 (save.punarg[ipPun][0] == 1 && strcmp(chTime,"LAST") != 0 ) )
3940 {
3941 fprintf( save.ipPnunit[ipPun],
3942 "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
3943 radius.Radius_mid_zone,
3944 radius.depth_mid_zone,
3945 wind.windv,
3946 wind.AccelTotalOutward,
3947 wind.AccelLine,
3948 wind.AccelCont ,
3949 wind.fmul ,
3950 wind.AccelGravity );
3951 }
3952 }
3953
3954 else if( strcmp(save.chSave[ipPun],"XATT") == 0 )
3955 {
3956 /* attenuated incident continuum */
3957 ASSERT( grid.lgOutputTypeOn[2] );
3958
3959 if( strcmp(chTime,"LAST") == 0 )
3960 {
3961 if( grid.lgGrid )
3962 saveFITSfile( save.ipPnunit[ipPun], 2 );
3963 else
3964 {
3965 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
3966 cdEXIT(EXIT_FAILURE);
3967 }
3968 }
3969 }
3970 else if( strcmp(save.chSave[ipPun],"XRFI") == 0 )
3971 {
3972 /* reflected incident continuum */
3973 ASSERT( grid.lgOutputTypeOn[3] );
3974
3975 if( strcmp(chTime,"LAST") == 0 )
3976 {
3977 if( grid.lgGrid )
3978 saveFITSfile( save.ipPnunit[ipPun], 3 );
3979 else
3980 {
3981 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
3982 cdEXIT(EXIT_FAILURE);
3983 }
3984 }
3985 }
3986 else if( strcmp(save.chSave[ipPun],"XINC") == 0 )
3987 {
3988 /* incident continuum */
3989 ASSERT( grid.lgOutputTypeOn[1] );
3990
3991 if( strcmp(chTime,"LAST") == 0 )
3992 {
3993 if( grid.lgGrid )
3994 saveFITSfile( save.ipPnunit[ipPun], 1 );
3995 else
3996 {
3997 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
3998 cdEXIT(EXIT_FAILURE);
3999 }
4000 }
4001 }
4002 else if( strcmp(save.chSave[ipPun],"XDFR") == 0 )
4003 {
4004 /* reflected diffuse continuous emission */
4005 ASSERT( grid.lgOutputTypeOn[5] );
4006
4007 if( strcmp(chTime,"LAST") == 0 )
4008 {
4009 if( grid.lgGrid )
4010 saveFITSfile( save.ipPnunit[ipPun], 5 );
4011 else
4012 {
4013 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4014 cdEXIT(EXIT_FAILURE);
4015 }
4016 }
4017 }
4018 else if( strcmp(save.chSave[ipPun],"XDFO") == 0 )
4019 {
4020 /* diffuse continuous emission outward */
4021 ASSERT( grid.lgOutputTypeOn[4] );
4022
4023 if( strcmp(chTime,"LAST") == 0 )
4024 {
4025 if( grid.lgGrid )
4026 saveFITSfile( save.ipPnunit[ipPun], 4 );
4027 else
4028 {
4029 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4030 cdEXIT(EXIT_FAILURE);
4031 }
4032 }
4033 }
4034 else if( strcmp(save.chSave[ipPun],"XLNR") == 0 )
4035 {
4036 /* reflected lines */
4037 ASSERT( grid.lgOutputTypeOn[7] );
4038
4039 if( strcmp(chTime,"LAST") == 0 )
4040 {
4041 if( grid.lgGrid )
4042 saveFITSfile( save.ipPnunit[ipPun], 7 );
4043 else
4044 {
4045 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4046 cdEXIT(EXIT_FAILURE);
4047 }
4048 }
4049 }
4050 else if( strcmp(save.chSave[ipPun],"XLNO") == 0 )
4051 {
4052 /* outward lines */
4053 ASSERT( grid.lgOutputTypeOn[6] );
4054
4055 if( strcmp(chTime,"LAST") == 0 )
4056 {
4057 if( grid.lgGrid )
4058 saveFITSfile( save.ipPnunit[ipPun], 6 );
4059 else
4060 {
4061 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4062 cdEXIT(EXIT_FAILURE);
4063 }
4064 }
4065 }
4066 else if( strcmp(save.chSave[ipPun],"XREF") == 0 )
4067 {
4068 /* total reflected, lines and continuum */
4069 ASSERT( grid.lgOutputTypeOn[9] );
4070
4071 if( strcmp(chTime,"LAST") == 0 )
4072 {
4073 if( grid.lgGrid )
4074 saveFITSfile( save.ipPnunit[ipPun], 9 );
4075 else
4076 {
4077 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4078 cdEXIT(EXIT_FAILURE);
4079 }
4080 }
4081 }
4082 else if( strcmp(save.chSave[ipPun],"XTOT") == 0 )
4083 {
4084 /* total spectrum, reflected plus transmitted */
4085 ASSERT( grid.lgOutputTypeOn[0] );
4086
4087 if( strcmp(chTime,"LAST") == 0 )
4088 {
4089 if( grid.lgGrid )
4090 saveFITSfile( save.ipPnunit[ipPun], 0 );
4091 else
4092 {
4093 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4094 cdEXIT(EXIT_FAILURE);
4095 }
4096 }
4097 }
4098 else if( strcmp(save.chSave[ipPun],"XTRN") == 0 )
4099 {
4100 /* total outward, lines and continuum */
4101 ASSERT( grid.lgOutputTypeOn[8] );
4102
4103 if( strcmp(chTime,"LAST") == 0 )
4104 {
4105 if( grid.lgGrid )
4106 saveFITSfile( save.ipPnunit[ipPun], 8 );
4107 else
4108 {
4109 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4110 cdEXIT(EXIT_FAILURE);
4111 }
4112 }
4113 }
4114 else if( strcmp(save.chSave[ipPun],"XSPM") == 0 )
4115 {
4116 /* exp(-tau) to the illuminated face */
4117 ASSERT( grid.lgOutputTypeOn[10] );
4118
4119 if( strcmp(chTime,"LAST") == 0 )
4120 {
4121 if( grid.lgGrid )
4122 saveFITSfile( save.ipPnunit[ipPun], 10 );
4123 else
4124 {
4125 fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4126 cdEXIT(EXIT_FAILURE);
4127 }
4128 }
4129 }
4130 // termination of second set of nested if's
4131 // error if we have not matched key
4132 /* there are a few "save" commands that are handled elsewhere
4133 * save dr is an example. These will have lgRealSave set false */
4134 // lgNoHitFirstBranch says did not find in previous nest of if's
4135 else if( save.lgRealSave[ipPun] && lgNoHitFirstBranch )
4136 {
4137 /* this is insanity, internal flag set in ParseSave not seen here */
4138 fprintf( ioQQQ, " PROBLEM DISASTER SaveDo does not recognize flag %4.4s set by ParseSave. This is impossible.\n",
4139 save.chSave[ipPun] );
4140 TotalInsanity();
4141 }
4142
4143 /* print special hash string to separate out various iterations
4144 * chTime is LAST on last iteration
4145 * save.lgHashEndIter flag is true by default, set false
4146 * with "no hash" keyword on save command
4147 * save.lg_separate_iterations is true by default, set false
4148 * when save time dependent calcs since do not want special
4149 * character between time steps
4150 * grid.lgGrid is only true when doing a grid of calculations */
4151 if( strcmp(chTime,"LAST") == 0 &&
4152 !(iterations.lgLastIt && !grid.lgGrid ) &&
4153 save.lgHashEndIter[ipPun] &&
4154 save.lg_separate_iterations[ipPun] &&
4155 !save.lgFITS[ipPun] )
4156 {
4157 if( dynamics.lgTimeDependentStatic && strcmp( save.chHashString , "TIME_DEP" )==0 )
4158 {
4159 fprintf( save.ipPnunit[ipPun], "\"time=%f\n",
4160 dynamics.time_elapsed );
4161 }
4162 else
4163 {
4164 fprintf( save.ipPnunit[ipPun], "%s",
4165 save.chHashString );
4166 if( grid.lgGrid && ( iterations.lgLastIt || lgAbort ) )
4167 fprintf( save.ipPnunit[ipPun], " GRID_DELIMIT -- grid%09ld",
4168 optimize.nOptimiz );
4169 fprintf( save.ipPnunit[ipPun], "\n" );
4170 }
4171 }
4172 if( save.lgFLUSH )
4173 fflush( save.ipPnunit[ipPun] );
4174 }
4175 }
4176 return;
4177}
4178
4179/*SaveLineIntensity produce the 'save lines intensity' output */
4180STATIC void SaveLineIntensity(FILE * ioPUN, long int ipPun , realnum Threshold )
4181{
4182 long int i;
4183
4184 DEBUG_ENTRY( "SaveLineIntensity()" );
4185
4186 /* used to save out all the emission line intensities
4187 * first initialize the line image reader */
4188
4189 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4190 input.echo(ioPUN);
4191
4192 /* now print any cautions or warnings */
4193 cdWarnings( ioPUN);
4194 cdCautions( ioPUN);
4195 fprintf( ioPUN, "zone=%5ld\n", nzone );
4196 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4197 fprintf( ioPUN, "begin emission lines\n" );
4198
4199 /* only save non-zero intensities */
4200 SaveResults1Line(ioPUN," ",0,0.,"Start");
4201
4202 // check whether intrinsic or emergent line emissivity
4203 bool lgEmergent = false;
4204 if( save.punarg[ipPun][0] > 0 )
4205 lgEmergent = true;
4206
4207 for( i=0; i < LineSave.nsum; i++ )
4208 {
4209 // Threshold is zero by default on save line intensity,
4210 // all option sets to negative number so that we report all lines
4211 if( LineSv[i].SumLine[lgEmergent] > Threshold )
4212 {
4213 SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
4214 LineSv[i].SumLine[save.lgEmergent[ipPun]], "Line ");
4215 }
4216 }
4217
4218 SaveResults1Line(ioPUN," ",0,0.,"Flush");
4219
4220 fprintf( ioPUN, " \n" );
4221 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4222
4223 return;
4224}
4225
4226/* lgSaveOpticalDepths true says save optical depths */
4227static bool lgPopsFirstCall , lgSaveOpticalDepths;
4228
4229/*SaveLineStuff save optical depths or source functions for all transferred lines */
4230STATIC void SaveLineStuff(
4231 FILE * ioPUN,
4232 const char *chJob ,
4233 realnum xLimit )
4234{
4235
4236 long int nelem , ipLo , ipHi , i , ipISO;
4237 long index = 0;
4238 static bool lgFirst=true;
4239
4240 DEBUG_ENTRY( "SaveLineStuff()" );
4241
4242 /*find out which job this is and set a flag to use later */
4243 if( strcmp( &*chJob , "optical" ) == 0 )
4244 {
4245 /* save line optical depths */
4246 lgSaveOpticalDepths = true;
4247 lgPopsFirstCall = false;
4248 }
4249 else if( strcmp( &*chJob , "populat" ) == 0 )
4250 {
4251 lgSaveOpticalDepths = false;
4252 /* level population information */
4253 if( lgFirst )
4254 {
4255 lgPopsFirstCall = true;
4256 fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
4257 lgFirst = false;
4258 }
4259 else
4260 {
4261 lgPopsFirstCall = false;
4262 }
4263 }
4264 else
4265 {
4266 fprintf( ioQQQ, " insane job in SaveLineStuff =%s\n",
4267 &*chJob );
4268 cdEXIT(EXIT_FAILURE);
4269 }
4270
4271 /* loop over all lines, calling put1Line to create info (routine located below) */
4272 /* hydrogen like lines */
4273 /* >>chng 02 may 16, had been explicit H and He-like loops */
4274 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
4275 {
4276 for( nelem=ipISO; nelem < LIMELM; nelem++ )
4277 {
4278 if( dense.lgElmtOn[nelem] )
4279 {
4280 /* 06 aug 28, from numLevels_max to _local. */
4281 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
4282 {
4283 for( ipLo=0; ipLo <ipHi; ipLo++ )
4284 {
4285 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
4286 continue;
4287
4288 ++index;
4289 Save1Line( iso_sp[ipISO][nelem].trans(ipHi,ipLo), ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
4290 }
4291 }
4292 /* also do extra Lyman lines if optical depths are to be done,
4293 * these are line that are included only for absorption, not in the
4294 * model atoms */
4295 if( lgSaveOpticalDepths )
4296 {
4297 /* >>chng 02 aug 23, for he-like, had starting on much too high a level since
4298 * index was number of levels - caught by Adrian Turner */
4299 /* now output extra line lines, starting one above those already done above */
4300 /*for( ipHi=iso_sp[ipISO][nelem].numLevels_max; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/
4301 /* 06 aug 28, from numLevels_max to _local. */
4302 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
4303 {
4304 ++index;
4305 Save1Line( ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
4306 }
4307 }
4308 }
4309 }
4310 }
4311
4312 /* index from 1 to skip over dummy line */
4313 for( i=1; i < nLevel1; i++ )
4314 {
4315 ++index;
4316 Save1Line( TauLines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
4317 }
4318
4319 for( i=0; i < nWindLine; i++ )
4320 {
4321 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
4322 {
4323 ++index;
4324 Save1Line( TauLine2[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
4325 }
4326 }
4327
4328 for( i=0; i < nUTA; i++ )
4329 {
4330 ++index;
4331 Save1Line( UTALines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
4332 }
4333
4334
4335 /* do optical depths of FeII lines, if large atom is turned on */
4336 FeIIPunchLineStuff( ioPUN , xLimit , index);
4337
4338 /* save optical depths of H2 lines */
4339 h2.H2_PunchLineStuff( ioPUN , xLimit , index);
4340
4341 /*fprintf(ioPUN, "##################################\n"); */
4342 fprintf( ioPUN , "%s\n",save.chHashString );
4343 return;
4344}
4345
4346/*Save1Line called by SaveLineStuff to produce output for one line */
4347void Save1Line( const TransitionProxy& t , FILE * ioPUN , realnum xLimit , long index, realnum DopplerWidth )
4348{
4349
4350 if( lgSaveOpticalDepths )
4351 {
4352 /* optical depths, no special first time, only print them */
4353 if( t.Emis().TauIn() >= xLimit )
4354 {
4355 /* label like "C 4" or "H 1" */
4356 fprintf( ioPUN , "%2.2s%2.2s\t",
4357 elementnames.chElementSym[(*t.Hi()).nelem()-1] ,
4358 elementnames.chIonStage[(*t.Hi()).IonStg()-1] );
4359
4360 /* print wavelengths, either line in main printout labels,
4361 * or in various units in exponential notation - prt_wl is in prt.c */
4362 if( strcmp( save.chConPunEnr[save.ipConPun], "labl" )== 0 )
4363 {
4364 prt_wl( ioPUN , t.WLAng() );
4365 }
4366 else
4367 {
4368 /* this converts energy in Rydbergs into any of the other units */
4369 fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t.EnergyRyd())) );
4370 }
4371 /* print the optical depth */
4372 fprintf( ioPUN , "\t%.3f", t.Emis().TauIn() );
4373 /* damping constant */
4374 fprintf(ioPUN, "\t%.3e",
4375 t.Emis().dampXvel() / DopplerWidth );
4376 fprintf(ioPUN, "\n");
4377 }
4378 }
4379 else if( lgPopsFirstCall )
4380 {
4381 char chLbl[11];
4382 /* first call to line populations, print atomic parameters and indices */
4383 strcpy( chLbl, chLineLbl(t) );
4384 fprintf(ioPUN, "%li\t%s" , index , chLbl );
4385 /* stat weights */
4386 fprintf(ioPUN, "\t%.0f\t%.0f",
4387 (*t.Lo()).g() ,(*t.Hi()).g());
4388 /* energy difference, gf */
4389 fprintf(ioPUN, "\t%.2f\t%.3e",
4390 t.EnergyWN() ,t.Emis().gf());
4391 fprintf(ioPUN, "\n");
4392 }
4393 else
4394 {
4395 /* not first call, so do level populations and indices defined above */
4396 if( (*t.Hi()).Pop() > xLimit )
4397 {
4398 /* >>chng 05 may 08, add abundances, which for iso-seq species is
4399 * the density of the parent ion, for other lines, is unity.
4400 * had not been included so pops for iso seq were rel to parent ion.
4401 * caught by John Everett */
4402 /* multiplication by abundance no longer necessary since iso pops now denormalized */
4403 fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index, (*t.Lo()).Pop(), (*t.Hi()).Pop() );
4404 }
4405 }
4406}
4407
4408/*SaveNewContinuum produce the 'save new continuum' output */
4409STATIC void SaveNewContinuum(FILE * ioPUN )
4410{
4411 long int ipLo, ipHi,
4412 j ,
4413 ncells;
4414
4415 double wllo, wlhi;
4416
4417 double *energy,
4418 *cont_incid,
4419 *cont_atten,
4420 *diffuse_in,
4421 *diffuse_out,
4422 *emis_lines_out,
4423 *emis_lines_in;
4424
4425 /* get the low limit */
4426 wllo = rfield.anu[0];
4427 /* get high-energy limit */
4428 wlhi = rfield.anu[rfield.nflux-1];
4429 /* use native continuum mesh */
4430 ipLo = ipoint(wllo)-1;
4431 ipHi = ipoint(wlhi)-1;
4432 ncells = ipHi - ipLo + 1;
4433
4434 /* now allocate the space */
4435 energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4436 cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4437 cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4438 diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4439 diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4440 emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4441 emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4442 /*emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));
4443 emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));*/
4444
4445 /* now convert to rydbergs */
4446 for( j=0; j<ncells; ++j )
4447 {
4448 energy[j] = rfield.AnuOrg[j+ipLo] - rfield.widflx[j+ipLo]/2.;
4449 }
4450
4451 fixit();
4452 /* all of these should be rerouted to cdSPEC2, but units not right
4453 * need ergs multiplied for one, and continuum and lines may not be added correctly.
4454 * Goal is to abandon current cdSPEC and replace it with current cdSPEC2 */
4455
4456#if 1
4457 /* for cdSPEC the energy vector is the lower edge of the energy cell */
4458 /* get incident continuum */
4459 cdSPEC( 1 , ncells , cont_incid );
4460 /* get attenuated incident continuum */
4461 cdSPEC( 2 , ncells , cont_atten );
4462 /* get diffuse continuous emission, reflected */
4463 cdSPEC( 5 , ncells , diffuse_in );
4464 /* get continuous emission outward direction */
4465 cdSPEC( 4 , ncells , diffuse_out );
4467 /* get all outward lines */
4468 cdSPEC( 6 , ncells , emis_lines_out );
4469 /* get all reflected lines */
4470 cdSPEC( 7 , ncells , emis_lines_in );
4471#else
4472 cdSPEC2( 1, rfield.nflux, 0, rfield.nflux - 1, cont_incid );
4473 /* get attenuated incident continuum */
4474 cdSPEC2( 2, rfield.nflux, 0, rfield.nflux - 1, cont_atten );
4475 /* get diffuse continuous emission, reflected */
4476 cdSPEC2( 5, rfield.nflux, 0, rfield.nflux - 1, diffuse_in );
4477 /* get continuous emission outward direction */
4478 cdSPEC2( 4, rfield.nflux, 0, rfield.nflux - 1, diffuse_out );
4479 /* get all outward lines */
4480 cdSPEC2( 6, rfield.nflux, 0, rfield.nflux - 1, emis_lines_out );
4481 /* get all reflected lines */
4482 cdSPEC2( 7, rfield.nflux, 0, rfield.nflux - 1, emis_lines_in );
4483#endif
4484
4485 /* for this example we will do a wavelength range */
4486 for( j=0; j<ncells-1; ++j )
4487 {
4488 /* photon energy in appropriate energy or wavelength units */
4489 fprintf( ioPUN,"%.5e\t", AnuUnit((realnum)(energy[j]+rfield.widflx[j+ipLo]/2.) ) );
4490 fprintf( ioPUN,"%.3e\t", cont_incid[j] );
4491 fprintf( ioPUN,"%.3e\t", cont_atten[j] );
4492 fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
4493 fprintf( ioPUN,"%.3e",
4494 emis_lines_out[j]+emis_lines_in[j]/*+emis_lines_pump_out[j]+emis_lines_pump_in[j]*/ );
4495 fprintf( ioPUN,"\n" );
4496 }
4497
4498 free(energy);
4499 free(cont_incid);
4500 free(diffuse_in);
4501 free(diffuse_out);
4502 free(cont_atten);
4503 free(emis_lines_out);
4504 free(emis_lines_in);
4505 /*free(emis_lines_pump_out);
4506 free(emis_lines_pump_in );*/
4507}
4508
4509/* save AGN3 hemiss, for Chapter 4, routine is below */
4510STATIC void AGN_Hemis(FILE *ioPUN )
4511{
4512 const int NTE = 4;
4513 realnum te[NTE] = {5000., 10000., 15000., 20000.};
4514 realnum *agn_continuum[NTE];
4515 double TempSave = phycon.te;
4516 long i , j;
4517
4518 DEBUG_ENTRY( "AGN_Hemis()" );
4519
4520 /* make table of continuous emission at various temperatuers */
4521 /* first allocate space */
4522 for( i=0;i<NTE; ++i)
4523 {
4524 agn_continuum[i] = (realnum *)MALLOC((unsigned)rfield.nflux*sizeof(realnum) );
4525
4526 /* set the next temperature */
4527 /* recompute everything at this new temp */
4528 TempChange(te[i] , true);
4529 /* converge the pressure-temperature-ionization solution for this zone */
4530 ConvPresTempEdenIoniz();
4531
4532 /* now get the thermal emission */
4533 RT_diffuse();
4534 for(j=0;j<rfield.nflux; ++j )
4535 {
4536 agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(realnum)dense.eden/
4537 (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] );
4538 }
4539 }
4540
4541 /* print title for line */
4542 fprintf(ioPUN,"wl");
4543 for( i=0;i<NTE; ++i)
4544 {
4545 fprintf(ioPUN,"\tT=%.0f",te[i]);
4546 }
4547 fprintf( ioPUN , "\tcont\n");
4548
4549 /* not print all n temperatures across a line */
4550 for(j=0;j<rfield.nflux; ++j )
4551 {
4552 fprintf( ioPUN , "%12.5e",
4553 AnuUnit(rfield.AnuOrg[j]) );
4554 /* loop over the temperatures, and for each, calculate a continuum */
4555 for( i=0;i<NTE; ++i)
4556 {
4557 fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
4558 }
4559 /* cont label and end of line*/
4560 fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j]);
4561 }
4562
4563 /* now free the continua */
4564 for( i=0;i<NTE; ++i)
4565 {
4566 free( agn_continuum[i] );
4567 }
4568
4569 /* Restore temperature stored before this routine was called */
4570 /* and force update */
4571 TempChange(TempSave , true);
4572
4573 fprintf( ioQQQ, "AGN_Hemis - result of save AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
4574 cdEXIT(EXIT_FAILURE);
4575}
4576
4577/*SaveResults save results from save results command */
4578/*SaveResults1Line do single line of output for the save results and save line intensity commands */
4579STATIC void SaveResults(FILE* ioPUN)
4580{
4581 long int i , nelem , ion;
4582
4583 DEBUG_ENTRY( "SaveResults()" );
4584
4585 /* used to save out line intensities, optical depths,
4586 * and column densities */
4587
4588 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4589 input.echo(ioPUN);
4590
4591 /* first print any cautions or warnings */
4592 cdWarnings(ioPUN);
4593 cdCautions(ioPUN);
4594 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4595
4596 fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );
4597
4598 fprintf( ioPUN, "BEGIN EMISSION LINES\n" );
4599 SaveResults1Line(ioPUN," ",0,0.,"Start");
4600
4601 for( i=0; i < LineSave.nsum; i++ )
4602 {
4603 if( LineSv[i].SumLine[0] > 0. )
4604 {
4605 SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
4606 LineSv[i].SumLine[0],"Line ");
4607 }
4608 }
4609
4610 SaveResults1Line(ioPUN," ",0,0.,"Flush");
4611
4612 fprintf( ioPUN, " \n" );
4613
4614 fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" );
4615
4616 /* this dumps out the whole array,*/
4617 /* following loop relies on LIMELM being 30, assert it here in case
4618 * this is ever changed */
4619 ASSERT( LIMELM == 30 );
4620 /* this order of indices is to keep 30 as the fastest variable,
4621 * and the 32 (LIMELM+1) as the slower one */
4622 for( nelem=0; nelem<LIMELM; nelem++ )
4623 {
4624 for(ion=0; ion < nelem+1; ion++)
4625 {
4626 fprintf( ioPUN, " %10.3e", mean.xIonMean[0][nelem][ion][0] );
4627 /* throw line feed every 10 numbers */
4628 if( nelem==9|| nelem==19 || nelem==29 )
4629 {
4630 fprintf( ioPUN, "\n" );
4631 }
4632 }
4633 }
4634
4635 fprintf( ioPUN, "END OF RESULTS\n" );
4636 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4637 return;
4638}
4639
4640static const int LINEWIDTH = 6;
4641
4642/*SaveResults1Line do single line of output for the save results and save line intensity commands */
4643/* the number of emission lines across one line of printout */
4644STATIC void SaveResults1Line(
4645 FILE * ioPUN,
4646 /* 4 char + null string */
4647 const char *chLab,
4648 realnum wl,
4649 double xInten,
4650 /* null terminated string saying what to do */
4651 const char *chFunction)
4652{
4653
4654 long int i;
4655 static realnum wavelength[LINEWIDTH];
4656 static long ipLine;
4657 static double xIntenSave[LINEWIDTH];
4658 static char chLabSave[LINEWIDTH][5];
4659
4660 DEBUG_ENTRY( "SaveResults1Line()" );
4661
4662 /* if LineWidth is changed then change format in write too */
4663
4664 if( strcmp(chFunction,"Start") == 0 )
4665 {
4666 ipLine = 0;
4667 }
4668 else if( strcmp(chFunction,"Line ") == 0 )
4669 {
4670 /* save results in array so that they can be printed when done */
4671 wavelength[ipLine] = wl;
4672 strcpy( chLabSave[ipLine], chLab );
4673 xIntenSave[ipLine] = xInten;
4674
4675 /* now increment the counter and then check if we have filled the line,
4676 * and so should print it */
4677 ++ipLine;
4678 /* do print now if we are in column mode (one line per line) or if we have filled up
4679 * the line */
4680 if( ( strcmp(save.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
4681 {
4682 /* "array " is usual array 6 wide, "column" is one line per line */
4683 for( i=0; i < ipLine; i++ )
4684 {
4685 fprintf( ioPUN, " %4.4s ", chLabSave[i] );
4686 prt_wl( ioPUN, wavelength[i] );
4687 fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
4688 /* >>chng 02 apr 24, do not print type */
4689 /* single column for input into data base */
4690 if( strcmp(save.chPunRltType,"column") == 0 )
4691 fprintf( ioPUN, "\n" );
4692 }
4693 /* only put cr if we did not just put one */
4694 if( strcmp(save.chPunRltType,"array ") == 0 )
4695 fprintf( ioPUN, " \n" );
4696 ipLine = 0;
4697 }
4698 }
4699 else if( strcmp(chFunction,"Flush") == 0 )
4700 {
4701 if( ipLine > 0 )
4702 {
4703 /* this is an option to print many emission lines across an output line,
4704 * the array option, or a single column of numbers, the column option
4705 * that appears on the "save results" and "save intensity" commands
4706 */
4707 /* usual array 6 wide */
4708 for( i=0; i < ipLine; i++ )
4709 {
4710 fprintf( ioPUN, " %4.4s", chLabSave[i] );
4711 prt_wl( ioPUN, wavelength[i] );
4712 fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
4713 /* >>chng 02 apr 24, do not print type */
4714 /* single column for input into data base */
4715 if( strcmp(save.chPunRltType,"column") == 0 )
4716 fprintf( ioPUN, "\n" );
4717 }
4718 if( strcmp(save.chPunRltType,"array ") == 0 )
4719 fprintf( ioPUN, " \n" );
4720 }
4721 }
4722 else
4723 {
4724 fprintf( ioQQQ, " SaveResults1Line called with insane request=%5.5s\n",
4725 chFunction );
4726 cdEXIT(EXIT_FAILURE);
4727 }
4728 return;
4729}
4730
4731static const int NENR_GAUNT = 37;
4732static const int NTE_GAUNT = 21;
4733
4734/*SaveGaunts called by save gaunts command to output gaunt factors */
4735STATIC void SaveGaunts(FILE* ioPUN)
4736{
4737 long int i,
4738 charge,
4739 ite,
4740 j;
4741
4742 realnum ener[NENR_GAUNT],
4743 ste[NTE_GAUNT],
4744 z;
4745 double tempsave;
4746 double g[NENR_GAUNT][NTE_GAUNT], gfac;
4747
4748
4749 DEBUG_ENTRY( "SaveGaunts()" );
4750
4751 /* this routine is called from the PUNCH GAUNTS command
4752 * it drives the gaunt factor routine to save gaunts over full range */
4753 tempsave = phycon.te;
4754
4755 for( i=0; i < NTE_GAUNT; i++ )
4756 {
4757 ste[i] = 0.5f*i;
4758 }
4759
4760 for( i=0; i < NENR_GAUNT; i++ )
4761 {
4762 ener[i] = 0.5f*i - 8.f;
4763 }
4764
4765 for( charge = 1; charge<LIMELM; charge++ )
4766 {
4767 /* nuclear charge */
4768 z = (realnum)log10((double)charge);
4769
4770 /* energy is log of energy */
4771 for( ite=0; ite < NTE_GAUNT; ite++ )
4772 {
4773 phycon.alogte = ste[ite];
4774 phycon.te = pow(10.,phycon.alogte);
4775
4776 for( j=0; j < NENR_GAUNT; j++ )
4777 {
4778 gfac = cont_gaunt_calc( phycon.te, (double)charge, pow(10.,(double)ener[j]) );
4779 /*fprintf(ioQQQ, "z %.2e ener %.2e temp %.2e gfac %.3e \n",
4780 pow(10.,z), pow(10.,ener[j]), (double)phycon.te, gfac );*/
4781 g[j][ite] = gfac;
4782 }
4783 }
4784
4785 /* now save out the results */
4786 fprintf( ioPUN, " " );
4787 for( i=1; i <= NTE_GAUNT; i++ )
4788 {
4789 fprintf( ioPUN, "\t%6.1f", ste[i-1] );
4790 }
4791 fprintf( ioPUN, "\n" );
4792
4793 for( j=0; j < NENR_GAUNT; j++ )
4794 {
4795 fprintf( ioPUN, "\t%6.2f", ener[j] );
4796 for( ite=0; ite < NTE_GAUNT; ite++ )
4797 {
4798 fprintf( ioPUN, "\t%6.2f", log10(g[j][ite]) );
4799 }
4800 fprintf( ioPUN, "\n" );
4801 }
4802
4803 fprintf( ioPUN, " " );
4804 for( i=0; i < NTE_GAUNT; i++ )
4805 {
4806 fprintf( ioPUN, "\t%6.1f", ste[i] );
4807 }
4808 fprintf( ioPUN, "\n" );
4809
4810 /* now do the same thing as triplets */
4811 fprintf( ioPUN, " " );
4812 for( i=0; i < NTE_GAUNT; i++ )
4813 {
4814 fprintf( ioPUN, "\t%6.1f", ste[i] );
4815 }
4816 fprintf( ioPUN, "\n" );
4817
4818 for( i=0; i < NTE_GAUNT; i++ )
4819 {
4820 for( j=0; j < NENR_GAUNT; j++ )
4821 {
4822 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", ste[i], ener[j],
4823 log10(g[j][i]) );
4824 }
4825 }
4826
4827 fprintf( ioPUN, "Below is log(u), log(gamma**2), gff\n" );
4828 /* do the same thing as above, but this time print log(u) and log(gamma2) instead of temp and energy. */
4829 for( i=0; i < NTE_GAUNT; i++ )
4830 {
4831 for( j=0; j < NENR_GAUNT; j++ )
4832 {
4833 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", 2.*z + log10(TE1RYD) - ste[i] , log10(TE1RYD)+ener[j]-ste[i],
4834 log10(g[j][i]) );
4835 }
4836 }
4837 fprintf( ioPUN, "end of charge = %li\n", charge);
4838 fprintf( ioPUN, "****************************\n");
4839 }
4840
4841 phycon.te = tempsave;
4842 phycon.alogte = log10( phycon.te );
4843 return;
4844}
4845
4846void SaveGrid(FILE* pnunit, exit_type status)
4847{
4848 if( pnunit == NULL )
4849 return;
4850
4851 if( optimize.nOptimiz == 0 )
4852 {
4853 /* start of line gives abort and warning summary */
4854 fprintf( pnunit, "#Index\tFailure?\tWarnings?\tExit code\t#rank\t#seq" );
4855 /* print start of each variable command line */
4856 for( int i=0; i < grid.nintparm; i++ )
4857 {
4858 char chStr[10];
4859 strncpy( chStr, optimize.chVarFmt[i], 9 );
4860 /* make sure this small bit of string is terminated */
4861 chStr[9] = '\0';
4862 fprintf( pnunit, "\t%s", chStr );
4863 }
4864 fprintf( pnunit, "\tgrid parameter string\n" );
4865 }
4866 /* abort / warning summary for this sim */
4867 bool lgNoFailure = ( status == ES_SUCCESS || status == ES_WARNINGS );
4868 fprintf( pnunit, "%9.9ld\t%c\t%c\t%20s\t%ld\t%ld",
4869 optimize.nOptimiz,
4870 TorF(!lgNoFailure),
4871 TorF(warnings.lgWarngs),
4872 cpu.i().chExitStatus(status).c_str(),
4873 cpu.i().nRANK(),
4874 grid.seqNum );
4875 /* the grid parameters */
4876 char chGridParam[INPUT_LINE_LENGTH];
4877 char chStringHold[100];
4878 sprintf( chStringHold, "%f", grid.interpParameters[optimize.nOptimiz][0] );
4879 strcpy( chGridParam, chStringHold );
4880 for( int j=0; j < grid.nintparm; j++ )
4881 {
4882 if( j > 0 )
4883 {
4884 sprintf( chStringHold, ", %f", grid.interpParameters[optimize.nOptimiz][j] );
4885 strcat( chGridParam, chStringHold );
4886 }
4887 fprintf( pnunit, "\t%f", grid.interpParameters[optimize.nOptimiz][j] );
4888 }
4889 fprintf( pnunit, "\t%s\n", chGridParam );
4890}
abund
t_abund abund
Definition: abund.cpp:5
ChargTranPun
void ChargTranPun(FILE *ipPnunit, char *chSave)
Definition: atmdat_char_tran.cpp:1729
CHIANTI_Upsilon
double CHIANTI_Upsilon(long, long, long, long, double)
Definition: species2.cpp:690
b1
static realnum b1[83]
Definition: atmdat_3body.cpp:25
a1
static double a1[83]
Definition: atmdat_3body.cpp:27
ipT157
long ipT157
Definition: atmdat_readin.cpp:42
ipT146
long ipT146
Definition: atmdat_readin.cpp:52
ipT370
long ipT370
Definition: atmdat_readin.cpp:42
ipT63
long ipT63
Definition: atmdat_readin.cpp:53
ipT610
long ipT610
Definition: atmdat_readin.cpp:41
ipT1032
long ipT1032
Definition: atmdat_readin.cpp:54
FeIIPunPop
void FeIIPunPop(FILE *ioPUN, bool lgPunchRange, long int ipRangeLo, long int ipRangeHi, bool lgPunchDensity)
Definition: atom_feii.cpp:2189
FeIIPunDepart
void FeIIPunDepart(FILE *ioPUN, bool lgDoAll)
Definition: atom_feii.cpp:2029
FeIIPunchLineStuff
void FeIIPunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: atom_feii.cpp:2820
FeIIPunchColden
void FeIIPunchColden(FILE *ioPUN)
Definition: atom_feii.cpp:1647
nFeIIConBins
long int nFeIIConBins
Definition: atom_feii.cpp:119
FeIIPunchLevels
void FeIIPunchLevels(FILE *ioPUN)
Definition: atom_feii.cpp:1625
FeIIPunchOpticalDepth
void FeIIPunchOpticalDepth(FILE *ioPUN)
Definition: atom_feii.cpp:1673
FeII_Cont
realnum ** FeII_Cont
Definition: atom_feii.cpp:113
FeIISaveLines
void FeIISaveLines(FILE *ioPUN)
Definition: atom_feii.cpp:1709
FeII
t_FeII FeII
Definition: atomfeii.cpp:5
called
t_called called
Definition: called.cpp:5
nzone
long int nzone
Definition: cddefines.cpp:14
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
iteration
long int iteration
Definition: cddefines.cpp:16
plankf
double plankf(long int ip)
Definition: service.cpp:1707
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
is_odd
bool is_odd(int j)
Definition: cddefines.h:714
ipSILICON
const int ipSILICON
Definition: cddefines.h:318
STATIC
#define STATIC
Definition: cddefines.h:97
ipRecEsc
const int ipRecEsc
Definition: cddefines.h:279
nMatch
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:451
POW3
#define POW3
Definition: cddefines.h:936
ipIRON
const int ipIRON
Definition: cddefines.h:330
exit_type
exit_type
Definition: cddefines.h:115
ES_WARNINGS
@ ES_WARNINGS
Definition: cddefines.h:118
ES_SUCCESS
@ ES_SUCCESS
Definition: cddefines.h:116
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
ipALUMINIUM
const int ipALUMINIUM
Definition: cddefines.h:317
MIN2
#define MIN2
Definition: cddefines.h:761
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
LIMELM
const int LIMELM
Definition: cddefines.h:258
POW2
#define POW2
Definition: cddefines.h:929
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:961
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
ipMAGNESIUM
const int ipMAGNESIUM
Definition: cddefines.h:316
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
ipSODIUM
const int ipSODIUM
Definition: cddefines.h:315
TorF
char TorF(bool l)
Definition: cddefines.h:710
AnuUnit
double AnuUnit(realnum energy)
Definition: service.cpp:173
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
NISO
const int NISO
Definition: cddefines.h:261
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
ipRecNetEsc
const int ipRecNetEsc
Definition: cddefines.h:281
realnum
float realnum
Definition: cddefines.h:103
MAX2
#define MAX2
Definition: cddefines.h:782
ipSULPHUR
const int ipSULPHUR
Definition: cddefines.h:320
max
long max(int a, long b)
Definition: cddefines.h:775
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
ipRecRad
const int ipRecRad
Definition: cddefines.h:283
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipNEON
const int ipNEON
Definition: cddefines.h:314
ShowMe
void ShowMe(void)
Definition: service.cpp:181
fixit
void fixit(void)
Definition: service.cpp:991
cdColm
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition: cddrive.cpp:636
cdCautions
void cdCautions(FILE *ioOUT)
Definition: cddrive.cpp:226
cdWarnings
void cdWarnings(FILE *ioPNT)
Definition: cddrive.cpp:198
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1228
cdExecTime
double cdExecTime()
Definition: cddrive.cpp:481
cdSPEC2
void cdSPEC2(int Option, long int nEnergy, long int ipLoEnergy, long int ipHiEnergy, realnum ReturnedSpectrum[])
cdSPEC
void cdSPEC(int Option, long int nEnergy, double ReturnedSpectrum[])
Definition: cdspec.cpp:21
LineSvSortWL
LinSv * LineSvSortWL
Definition: cdinit.cpp:71
nWindLine
long nWindLine
Definition: cdinit.cpp:19
LineSv
LinSv * LineSv
Definition: cdinit.cpp:70
CollisionProxy::ColUL
realnum ColUL(const ColliderList &colls) const
Definition: collision.h:99
EmissionProxy::damp
realnum & damp() const
Definition: emission.h:563
EmissionProxy::TauCon
realnum & TauCon() const
Definition: emission.h:453
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:523
EmissionProxy::Pelec_esc
realnum & Pelec_esc() const
Definition: emission.h:533
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:423
EmissionProxy::pump
double & pump() const
Definition: emission.h:473
EmissionProxy::Pdest
realnum & Pdest() const
Definition: emission.h:543
EmissionProxy::TauTot
realnum & TauTot() const
Definition: emission.h:433
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
Energy::Ryd
double Ryd() const
Definition: energy.h:26
GrainVar::TotalEden
double TotalEden
Definition: grainvar.h:528
GrainVar::GraphiteEmission
vector< realnum > GraphiteEmission
Definition: grainvar.h:579
GrainVar::dstsc
vector< double > dstsc
Definition: grainvar.h:525
GrainVar::GrainEmission
vector< realnum > GrainEmission
Definition: grainvar.h:578
GrainVar::dstab
vector< double > dstab
Definition: grainvar.h:524
GrainVar::rate_h2_form_grains_used_total
double rate_h2_form_grains_used_total
Definition: grainvar.h:574
GrainVar::GasCoolColl
double GasCoolColl
Definition: grainvar.h:544
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
GrainVar::SilicateEmission
vector< realnum > SilicateEmission
Definition: grainvar.h:580
Singleton< t_yield >::Inst
static t_yield & Inst()
Definition: cddefines.h:175
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
TransitionProxy::Coll
CollisionProxy Coll() const
Definition: transition.h:424
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
TransitionProxy::EnergyRyd
double EnergyRyd() const
Definition: transition.h:83
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:450
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
diatomics::H2_PunchDo
void H2_PunchDo(FILE *io, char chJOB[], const char chTime[], long int ipPun)
Definition: mole_h2_io.cpp:1224
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:345
diatomics::H2_PunchLineStuff
void H2_PunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: mole_h2_io.cpp:1167
molecule
Definition: mole.h:132
molecule::index
int index
Definition: mole.h:169
molezone::den
double den
Definition: mole.h:358
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat.h:329
t_cpu_i::lgMPI_talk
bool lgMPI_talk() const
Definition: cpu.h:328
t_cpu_i::chExitStatus
const string & chExitStatus(exit_type s) const
Definition: cpu.h:333
t_cpu_i::big_endian
bool big_endian() const
Definition: cpu.h:287
t_cpu_i::nRANK
long nRANK() const
Definition: cpu.h:326
t_cpu::i
t_cpu_i & i()
Definition: cpu.h:347
t_ionbal::CollIonRate_Ground
double *** CollIonRate_Ground
Definition: ionbal.h:120
t_ionbal::RR_rate_coef_used
double ** RR_rate_coef_used
Definition: ionbal.h:212
t_ionbal::UTA_ionize_rate
double ** UTA_ionize_rate
Definition: ionbal.h:170
t_ionbal::CotaRate
realnum CotaRate[LIMELM]
Definition: ionbal.h:242
t_ionbal::PhotoRate_Shell
double **** PhotoRate_Shell
Definition: ionbal.h:111
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:198
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion)
Definition: ionbal.h:254
t_ionbal::DR_Badnell_rate_coef
double ** DR_Badnell_rate_coef
Definition: ionbal.h:205
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_iso_sp::RadRec_caseB
double RadRec_caseB
Definition: iso.h:513
t_iso_sp::st
qList st
Definition: iso.h:453
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:334
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
t_mole_global::list
MoleculeList list
Definition: mole.h:317
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
t_mole_local::findrk
double findrk(const char buf[]) const
Definition: mole_reactions.cpp:3886
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
t_yield::nelec_eject
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
t_yield::elec_eject_frac
realnum elec_eject_frac(long n, long i, long ns, long ne) const
Definition: yield.h:48
colden
t_colden colden
Definition: colden.cpp:5
ipCOL_elec
#define ipCOL_elec
Definition: colden.h:30
ipCOL_HTOT
#define ipCOL_HTOT
Definition: colden.h:12
ipCOL_H2s
#define ipCOL_H2s
Definition: colden.h:18
ipCOL_H2g
#define ipCOL_H2g
Definition: colden.h:16
ipCOL_Hp
#define ipCOL_Hp
Definition: colden.h:26
ipCOL_H0
#define ipCOL_H0
Definition: colden.h:22
ipCOL_H3p
#define ipCOL_H3p
Definition: colden.h:28
colliders
ColliderList colliders
Definition: collision.cpp:57
ipELECTRON
@ ipELECTRON
Definition: collision.h:9
GammaPrt
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
Definition: cont_gammas.cpp:253
cont_gaunt_calc
double cont_gaunt_calc(double temp, double z, double photon)
Definition: cont_gaunt.cpp:26
ipFineCont
long ipFineCont(double energy_ryd)
Definition: cont_ipoint.cpp:226
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:16
continuum
t_continuum continuum
Definition: continuum.cpp:5
conv
t_conv conv
Definition: conv.cpp:5
ConvPresTempEdenIoniz
int ConvPresTempEdenIoniz(void)
Definition: conv_pres_temp_eden_ioniz.cpp:23
NTYPES
@ NTYPES
Definition: conv.h:83
CoolSave
void CoolSave(FILE *io, char chJob[])
Definition: cool_save.cpp:20
cpu
static t_cpu cpu
Definition: cpu.h:355
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
dense
t_dense dense
Definition: dense.cpp:24
DoppVel
t_DoppVel DoppVel
Definition: doppvel.cpp:5
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:499
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:44
DynaSave
void DynaSave(FILE *ipPnunit, char chJob)
Definition: dynamics.cpp:2148
DynaPunchTimeDep
void DynaPunchTimeDep(FILE *ipPnunit, const char *chJob)
Definition: dynamics.cpp:2039
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
geometry
t_geometry geometry
Definition: geometry.cpp:5
gv
GrainVar gv
Definition: grainvar.cpp:5
grid
t_grid grid
Definition: grid.cpp:5
GridGatherInCloudy
void GridGatherInCloudy(void)
Definition: grid_xspec.cpp:167
NUM_OUTPUT_TYPES
const int NUM_OUTPUT_TYPES
Definition: grid.h:21
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
hcmap
t_hcmap hcmap
Definition: hcmap.cpp:21
map_do
void map_do(FILE *io, const char *chType)
Definition: hcmap.cpp:23
SaveHeat
void SaveHeat(FILE *io)
Definition: heat_save.cpp:22
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
AGN_He1_CS
void AGN_He1_CS(FILE *ioPun)
Definition: iso_solve.cpp:464
hmi
t_hmi hmi
Definition: hmi.cpp:5
HydroRecCool
double HydroRecCool(long int n, long int ipZ)
Definition: hydroreccool.cpp:10
hyperfine
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
input
t_input input
Definition: input.cpp:12
ion_recombAGN
void ion_recombAGN(FILE *io)
Definition: ion_recomb.cpp:218
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
ipH1s
const int ipH1s
Definition: iso.h:27
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:49
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
ipH2p
const int ipH2p
Definition: iso.h:29
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
ipH2s
const int ipH2s
Definition: iso.h:28
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
iterations
t_iterations iterations
Definition: iterations.cpp:5
LineSave
t_LineSave LineSave
Definition: lines.cpp:5
lines_service.h
magnetic
t_magnetic magnetic
Definition: magnetic.cpp:17
mean
t_mean mean
Definition: mean.cpp:17
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
mole
t_mole_local mole
Definition: mole.cpp:7
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
mole_punch
void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth)
Definition: mole_reactions.cpp:4221
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:814
null_mole
molecule * null_mole
Definition: mole_species.cpp:65
mole_dominant_rates
void mole_dominant_rates(const molecule *debug_species, FILE *ioOut)
Definition: mole_eval_balance.cpp:335
mole_priv.h
wavelength
static realnum * wavelength
Definition: monitor_results.cpp:70
lgCheckMonitors
bool lgCheckMonitors(FILE *ioMONITOR)
Definition: monitor_results.cpp:1603
monitor_results.h
opac
t_opac opac
Definition: opacity.cpp:5
optimize
t_optimize optimize
Definition: optimize.cpp:5
phycon
t_phycon phycon
Definition: phycon.cpp:6
BOLTZMANN
UNUSED const double BOLTZMANN
Definition: physconst.h:97
EN1RYD
UNUSED const double EN1RYD
Definition: physconst.h:179
EVRYD
UNUSED const double EVRYD
Definition: physconst.h:189
TE1RYD
UNUSED const double TE1RYD
Definition: physconst.h:183
RYDLAM
UNUSED const double RYDLAM
Definition: physconst.h:176
pressure
t_pressure pressure
Definition: pressure.cpp:5
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
prt_LineLabels
void prt_LineLabels(FILE *ioOUT, bool lgPrintAll)
Definition: prt.cpp:168
sprt_wl
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:25
PrtMeanIon
void PrtMeanIon(char chType, bool lgDensity, FILE *)
Definition: prt_meanion.cpp:11
PrtLinePres
void PrtLinePres(FILE *ioPRESSURE)
Definition: prt_linepres.cpp:16
PrtColumns
void PrtColumns(FILE *ioMEAN, const char *chType, long int ipPun)
Definition: prt_columns.cpp:14
lgFirst
static bool lgFirst
Definition: prt_linesum.cpp:16
ipLine
static long int * ipLine
Definition: prt_linesum.cpp:15
radius
t_radius radius
Definition: radius.cpp:5
rfield
t_rfield rfield
Definition: rfield.cpp:8
RT_diffuse
void RT_diffuse(void)
Definition: rt_diffuse.cpp:34
save
t_save save
Definition: save.cpp:5
MAX_HEADER_SIZE
static const long MAX_HEADER_SIZE
Definition: save.h:12
save_line
void save_line(FILE *ip, const char *chDo, bool lgEmergent)
Definition: save_line.cpp:102
VERSION_TRNCON
static const long VERSION_TRNCON
Definition: save.h:15
save_colden
void save_colden(FILE *ioPUN)
Definition: save_colden.cpp:92
SaveSpecial
void SaveSpecial(FILE *io, const char *chTime)
Definition: save_special.cpp:14
save_average
void save_average(long int ipPun)
Definition: save_average.cpp:210
SaveSpecies
void SaveSpecies(FILE *ioPUN, long int ipPun)
Definition: save_species.cpp:17
SaveLineData
NORETURN void SaveLineData(FILE *io)
Definition: save_linedata.cpp:25
Save_Line_RT
void Save_Line_RT(FILE *ip)
Definition: save_line.cpp:351
save_opacity
void save_opacity(FILE *io, long int np)
Definition: save_opacity.cpp:22
saveFITSfile
void saveFITSfile(FILE *io, int option)
Definition: save_fits.cpp:85
AGN_Hemis
STATIC void AGN_Hemis(FILE *ioPUN)
Definition: save_do.cpp:4510
wavelength_compare
int wavelength_compare(const void *a, const void *b)
Definition: save_do.cpp:75
lgSaveOpticalDepths
static bool lgSaveOpticalDepths
Definition: save_do.cpp:4227
Save1Line
void Save1Line(const TransitionProxy &t, FILE *ioPUN, realnum xLimit, long index, realnum DopplerWidth)
Definition: save_do.cpp:4347
lgPopsFirstCall
static bool lgPopsFirstCall
Definition: save_do.cpp:4227
LINEWIDTH
static const int LINEWIDTH
Definition: save_do.cpp:4640
SaveResults1Line
STATIC void SaveResults1Line(FILE *ioPUN, const char *chLab, realnum wl, double xInten, const char *chFunction)
Definition: save_do.cpp:4644
SaveGrid
void SaveGrid(FILE *pnunit, exit_type status)
Definition: save_do.cpp:4846
SaveFeII_cont
STATIC realnum SaveFeII_cont(long int ipCont, long ipFeII_Cont_type)
Definition: save_do.cpp:561
SaveLineStuff
STATIC void SaveLineStuff(FILE *ioPUN, const char *chJob, realnum xLimit)
Definition: save_do.cpp:4230
FindStrongestLineLabels
STATIC void FindStrongestLineLabels(void)
Definition: save_do.cpp:89
PrettyTransmission
realnum PrettyTransmission(long j, realnum transmission)
Definition: save_do.cpp:515
NTE_GAUNT
static const int NTE_GAUNT
Definition: save_do.cpp:4732
SaveDo
void SaveDo(const char *chTime)
Definition: save_do.cpp:573
SaveLineIntensity
STATIC void SaveLineIntensity(FILE *ioPUN, long int ipPun, realnum Threshold)
Definition: save_do.cpp:4180
chDummy
char * chDummy
Definition: save_do.cpp:558
SaveNewContinuum
STATIC void SaveNewContinuum(FILE *ioPUN)
Definition: save_do.cpp:4409
SaveGaunts
STATIC void SaveGaunts(FILE *ioPUN)
Definition: save_do.cpp:4735
SaveResults
STATIC void SaveResults(FILE *ioPUN)
Definition: save_do.cpp:4579
NENR_GAUNT
static const int NENR_GAUNT
Definition: save_do.cpp:4731
lgMustPrintHeader
static bool lgMustPrintHeader
Definition: save_line.cpp:287
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
g
static double * g
Definition: species2.cpp:28
struc
t_struc struc
Definition: struc.cpp:6
Wind::fmul
realnum fmul
Definition: wind.h:65
Wind::windv
realnum windv
Definition: wind.h:18
Wind::dvdr
realnum dvdr
Definition: wind.h:21
Wind::AccelTotalOutward
realnum AccelTotalOutward
Definition: wind.h:52
Wind::AccelLine
realnum AccelLine
Definition: wind.h:61
Wind::AccelCont
realnum AccelCont
Definition: wind.h:55
Wind::AccelGravity
realnum AccelGravity
Definition: wind.h:49
t_DoppVel::lgTurb_pressure
bool lgTurb_pressure
Definition: doppvel.h:33
t_DoppVel::TurbVel
realnum TurbVel
Definition: doppvel.h:12
t_FeII::nFeIILevel_malloc
long int nFeIILevel_malloc
Definition: atomfeii.h:193
t_FeII::FeIIAul
double FeIIAul[NFE2LEVN][NFE2LEVN]
Definition: atomfeii.h:213
t_FeII::FeIINRGs
double FeIINRGs[NFE2LEVN]
Definition: atomfeii.h:207
t_FeII::FeIISTWT
double FeIISTWT[NFE2LEVN]
Definition: atomfeii.h:210
t_FeII::FeIIColl
double FeIIColl[NFE2LEVN][NFE2LEVN][8]
Definition: atomfeii.h:216
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_Heavy::chShell
char chShell[7][3]
Definition: heavy.h:31
t_Heavy::Valence_IP_Ryd
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
t_LineSave::nsum
long int nsum
Definition: lines.h:62
t_LineSave::ipass
long int ipass
Definition: lines.h:75
t_called::lgTalk
bool lgTalk
Definition: called.h:12
t_colden::H0_ov_Tspin
double H0_ov_Tspin
Definition: colden.h:58
t_colden::C2Pops
realnum C2Pops[5]
Definition: colden.h:64
t_colden::H0_21cm_lower
double H0_21cm_lower
Definition: colden.h:96
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:38
t_colden::C1Pops
realnum C1Pops[3]
Definition: colden.h:76
t_colden::O1Pops
realnum O1Pops[3]
Definition: colden.h:80
t_colden::H0_21cm_upper
double H0_21cm_upper
Definition: colden.h:95
t_continuum::ResolutionScaleFactor
double ResolutionScaleFactor
Definition: continuum.h:90
t_continuum::mesh_md5sum
string mesh_md5sum
Definition: continuum.h:127
t_conv::lgConvTemp
bool lgConvTemp
Definition: conv.h:196
t_conv::hist_pres_density
vector< double > hist_pres_density
Definition: conv.h:295
t_conv::hist_pres_error
vector< double > hist_pres_error
Definition: conv.h:295
t_conv::hist_pres_current
vector< double > hist_pres_current
Definition: conv.h:295
t_conv::nPres2Ioniz
long int nPres2Ioniz
Definition: conv.h:152
t_conv::getCounterZone
long getCounterZone(const long type)
Definition: conv.h:327
t_conv::hist_temp_heat
vector< double > hist_temp_heat
Definition: conv.h:300
t_conv::lgConvPres
bool lgConvPres
Definition: conv.h:199
t_conv::hist_temp_cool
vector< double > hist_temp_cool
Definition: conv.h:300
t_conv::hist_temp_temp
vector< double > hist_temp_temp
Definition: conv.h:300
t_dense::pden
realnum pden
Definition: dense.h:98
t_dense::eden
double eden
Definition: dense.h:190
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:216
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
t_dense::EdenTrue
double EdenTrue
Definition: dense.h:221
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
t_dense::xMassDensity
realnum xMassDensity
Definition: dense.h:91
t_dynamics::Cool
double Cool()
Definition: dynamics.cpp:2187
t_dynamics::lgTimeDependentStatic
bool lgTimeDependentStatic
Definition: dynamics.h:96
t_dynamics::Source
double ** Source
Definition: dynamics.h:74
t_dynamics::Rate
double Rate
Definition: dynamics.h:71
t_dynamics::Heat
double Heat()
Definition: dynamics.cpp:2173
t_dynamics::time_elapsed
double time_elapsed
Definition: dynamics.h:99
t_dynamics::lgStatic_completed
bool lgStatic_completed
Definition: dynamics.h:105
t_elementnames::chIonStage
char chIonStage[LIMELM+1][CHARS_ION_STAGE]
Definition: elementnames.h:29
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_geometry::lgSphere
bool lgSphere
Definition: geometry.h:24
t_geometry::FillFac
realnum FillFac
Definition: geometry.h:19
t_geometry::covgeo
realnum covgeo
Definition: geometry.h:35
t_grid::lgGrid
bool lgGrid
Definition: grid.h:40
t_grid::seqNum
long seqNum
Definition: grid.h:62
t_grid::nintparm
long nintparm
Definition: grid.h:47
t_grid::lgOutputTypeOn
bool lgOutputTypeOn[NUM_OUTPUT_TYPES]
Definition: grid.h:56
t_grid::interpParameters
realnum ** interpParameters
Definition: grid.h:31
t_hcmap::MapZone
long int MapZone
Definition: hcmap.h:20
t_hcmap::lgMapDone
bool lgMapDone
Definition: hcmap.h:36
t_hcmap::lgMapBeingDone
bool lgMapBeingDone
Definition: hcmap.h:33
t_hmi::HeatH2Dexc_used
double HeatH2Dexc_used
Definition: hmi.h:137
t_hmi::UV_Cont_rel2_Draine_DB96_depth
realnum UV_Cont_rel2_Draine_DB96_depth
Definition: hmi.h:74
t_hmi::HeatH2Dish_used
double HeatH2Dish_used
Definition: hmi.h:129
t_hmi::H2_Solomon_dissoc_rate_used_H2g
double H2_Solomon_dissoc_rate_used_H2g
Definition: hmi.h:92
t_hmi::H2_total
double H2_total
Definition: hmi.h:16
t_hmi::H2_photodissoc_used_H2g
double H2_photodissoc_used_H2g
Definition: hmi.h:108
t_hmi::UV_Cont_rel2_Habing_TH85_depth
realnum UV_Cont_rel2_Habing_TH85_depth
Definition: hmi.h:64
t_hyperfine::Tspin21cm
double Tspin21cm
Definition: hyperfine.h:47
t_input::echo
void echo(FILE *ipOUT)
Definition: input.cpp:89
t_iterations::lgLastIt
bool lgLastIt
Definition: iterations.h:36
t_magnetic::EnthalpyDensity
double EnthalpyDensity
Definition: magnetic.h:30
t_magnetic::pressure
double pressure
Definition: magnetic.h:33
t_mean::xIonMean
multi_arr< double, 4 > xIonMean
Definition: mean.h:14
t_opac::TauAbsFace
realnum * TauAbsFace
Definition: opacity.h:91
t_opac::ipElement
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:269
t_opac::opacity_sct
double * opacity_sct
Definition: opacity.h:98
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:95
t_opac::TauAbsGeo
realnum ** TauAbsGeo
Definition: opacity.h:82
t_opac::telec
realnum telec
Definition: opacity.h:175
t_opac::TauScatFace
realnum * TauScatFace
Definition: opacity.h:92
t_opac::albedo
double * albedo
Definition: opacity.h:104
t_optimize::nOptimiz
long int nOptimiz
Definition: optimize.h:246
t_optimize::chVarFmt
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:263
t_phycon::te_eV
double te_eV
Definition: phycon.h:14
t_phycon::EnergyExcitation
double EnergyExcitation
Definition: phycon.h:37
t_phycon::te03
double te03
Definition: phycon.h:59
t_phycon::EnergyBinding
double EnergyBinding
Definition: phycon.h:44
t_phycon::te
double te
Definition: phycon.h:11
t_phycon::EnthalpyDensity
double EnthalpyDensity
Definition: phycon.h:40
t_phycon::telogn
double telogn[7]
Definition: phycon.h:76
t_phycon::TeProp
double TeProp
Definition: phycon.h:91
t_phycon::alogte
double alogte
Definition: phycon.h:82
t_phycon::te70
double te70
Definition: phycon.h:51
t_phycon::TeInit
double TeInit
Definition: phycon.h:89
t_phycon::EnergyIonization
double EnergyIonization
Definition: phycon.h:31
t_phycon::te10
double te10
Definition: phycon.h:55
t_pressure::PresTotlError
double PresTotlError
Definition: pressure.h:87
t_pressure::PresGasCurr
double PresGasCurr
Definition: pressure.h:89
t_pressure::IntegRhoGravity
double IntegRhoGravity
Definition: pressure.h:123
t_pressure::PresRamCurr
double PresRamCurr
Definition: pressure.h:79
t_pressure::PresTurbCurr
double PresTurbCurr
Definition: pressure.h:82
t_pressure::pinzon_PresIntegElecThin
realnum pinzon_PresIntegElecThin
Definition: pressure.h:116
t_pressure::PresInteg
realnum PresInteg
Definition: pressure.h:109
t_pressure::pinzon
realnum pinzon
Definition: pressure.h:110
t_pressure::PresIntegElecThin
realnum PresIntegElecThin
Definition: pressure.h:115
t_pressure::pres_radiation_lines_curr
double pres_radiation_lines_curr
Definition: pressure.h:101
t_pressure::PresTotlInit
double PresTotlInit
Definition: pressure.h:92
t_pressure::PresTotlCurr
double PresTotlCurr
Definition: pressure.h:86
t_radius::drad
double drad
Definition: radius.h:31
t_radius::r1r0sq
double r1r0sq
Definition: radius.h:49
t_radius::depth_mid_zone
double depth_mid_zone
Definition: radius.h:41
t_radius::drad_x_fillfac
double drad_x_fillfac
Definition: radius.h:71
t_radius::Conv2PrtInten
double Conv2PrtInten
Definition: radius.h:147
t_radius::Radius_mid_zone
double Radius_mid_zone
Definition: radius.h:28
t_radius::Radius
double Radius
Definition: radius.h:25
t_radius::depth
double depth
Definition: radius.h:38
t_rfield::extin_mag_V_extended
double extin_mag_V_extended
Definition: rfield.h:281
t_rfield::comup
double * comup
Definition: rfield.h:255
t_rfield::ConEmitReflec
realnum ** ConEmitReflec
Definition: rfield.h:155
t_rfield::fine_opt_depth
realnum * fine_opt_depth
Definition: rfield.h:410
t_rfield::ConEmitLocal
realnum ** ConEmitLocal
Definition: rfield.h:149
t_rfield::outlin_noplot
realnum * outlin_noplot
Definition: rfield.h:200
t_rfield::chContLabel
char ** chContLabel
Definition: rfield.h:223
t_rfield::emm
realnum emm
Definition: rfield.h:49
t_rfield::nupper
long int nupper
Definition: rfield.h:46
t_rfield::getCoarseTransCoef
const realnum * getCoarseTransCoef()
Definition: rfield.cpp:10
t_rfield::fine_anu
realnum * fine_anu
Definition: rfield.h:412
t_rfield::flux
realnum ** flux
Definition: rfield.h:86
t_rfield::flux_total_incident
realnum ** flux_total_incident
Definition: rfield.h:209
t_rfield::ipV_filter
long int ipV_filter
Definition: rfield.h:259
t_rfield::TotDiff2Pht
realnum * TotDiff2Pht
Definition: rfield.h:187
t_rfield::ipMaxBolt
long int ipMaxBolt
Definition: rfield.h:249
t_rfield::extin_mag_V_point
double extin_mag_V_point
Definition: rfield.h:277
t_rfield::reflin
realnum ** reflin
Definition: rfield.h:206
t_rfield::otscon
realnum * otscon
Definition: rfield.h:195
t_rfield::anu2
realnum * anu2
Definition: rfield.h:79
t_rfield::chLineLabel
char ** chLineLabel
Definition: rfield.h:220
t_rfield::ConEmitOut
realnum ** ConEmitOut
Definition: rfield.h:161
t_rfield::otslin
realnum * otslin
Definition: rfield.h:193
t_rfield::line_count
long int * line_count
Definition: rfield.h:68
t_rfield::AnuOrg
double * AnuOrg
Definition: rfield.h:62
t_rfield::nfine
long nfine
Definition: rfield.h:402
t_rfield::egamry
realnum egamry
Definition: rfield.h:52
t_rfield::outlin
realnum ** outlin
Definition: rfield.h:199
t_rfield::comdn
double * comdn
Definition: rfield.h:256
t_rfield::nflux
long int nflux
Definition: rfield.h:43
t_rfield::widflx
realnum * widflx
Definition: rfield.h:65
t_rfield::trans_coef_total
realnum * trans_coef_total
Definition: rfield.h:390
t_rfield::DiffuseLineEmission
realnum * DiffuseLineEmission
Definition: rfield.h:203
t_rfield::anu
double * anu
Definition: rfield.h:58
t_rfield::ConSourceFcnLocal
realnum ** ConSourceFcnLocal
Definition: rfield.h:152
t_rfield::fine_opac_zone
realnum * fine_opac_zone
Definition: rfield.h:408
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:164
t_rfield::ConRefIncid
realnum ** ConRefIncid
Definition: rfield.h:167
t_save::lgFLUSH
bool lgFLUSH
Definition: save.h:298
t_save::chHeader
char chHeader[LIMPUN][MAX_HEADER_SIZE]
Definition: save.h:231
t_save::chLineListLabel
vector< char * > chLineListLabel[LIMPUN]
Definition: save.h:178
t_save::chPunRltType
char chPunRltType[7]
Definition: save.h:319
t_save::ResolutionAbs
realnum ResolutionAbs
Definition: save.h:366
t_save::chConPunEnr
const char * chConPunEnr[LIMPUN]
Definition: save.h:284
t_save::ipConPun
long int ipConPun
Definition: save.h:287
t_save::lg_separate_iterations
bool lg_separate_iterations[LIMPUN]
Definition: save.h:242
t_save::lgEmergent
bool lgEmergent[LIMPUN]
Definition: save.h:216
t_save::nSaveEveryZone
long int nSaveEveryZone[LIMPUN]
Definition: save.h:262
t_save::LinEvery
long int LinEvery
Definition: save.h:350
t_save::lgPunHeader
bool lgPunHeader[LIMPUN]
Definition: save.h:246
t_save::punarg
realnum punarg[LIMPUN][3]
Definition: save.h:254
t_save::whichDiatomToPrint
diatomics * whichDiatomToPrint[LIMPUN]
Definition: save.h:226
t_save::ncSaveSkip
long int ncSaveSkip
Definition: save.h:354
t_save::chNONSENSE
const char * chNONSENSE
Definition: save.h:233
t_save::chSave
char chSave[LIMPUN][5]
Definition: save.h:225
t_save::nsave
long int nsave
Definition: save.h:222
t_save::lgCumulative
bool lgCumulative[LIMPUN]
Definition: save.h:219
t_save::lgHashEndIter
bool lgHashEndIter[LIMPUN]
Definition: save.h:291
t_save::lgLinEvery
bool lgLinEvery
Definition: save.h:351
t_save::optname
string optname[LIMPUN]
Definition: save.h:257
t_save::chHashString
char chHashString[INPUT_LINE_LENGTH]
Definition: save.h:295
t_save::nLineList
long nLineList[LIMPUN]
Definition: save.h:176
t_save::chSpeciesDominantRates
char chSpeciesDominantRates[LIMPUN][CHARS_SPECIES]
Definition: save.h:368
t_save::lgPunLstIter
bool lgPunLstIter[LIMPUN]
Definition: save.h:271
t_save::lgRealSave
bool lgRealSave[LIMPUN]
Definition: save.h:213
t_save::chSaveArgs
char chSaveArgs[LIMPUN][5]
Definition: save.h:265
t_save::lgSaveEveryZone
bool lgSaveEveryZone[LIMPUN]
Definition: save.h:261
t_save::ipPnunit
FILE * ipPnunit[LIMPUN]
Definition: save.h:197
t_save::emisfreq
Energy emisfreq[LIMPUN]
Definition: save.h:371
t_save::ipEmisFreq
long ipEmisFreq[LIMPUN]
Definition: save.h:372
t_save::lgLineListRatio
bool lgLineListRatio[LIMPUN]
Definition: save.h:182
t_save::Resolution
realnum Resolution
Definition: save.h:363
t_save::wlLineList
vector< realnum > wlLineList[LIMPUN]
Definition: save.h:180
t_save::lgFITS
bool lgFITS[LIMPUN]
Definition: save.h:274
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
t_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:53
t_species::numLevels_max
long numLevels_max
Definition: cddefines.h:1239
t_species::chLabel
char * chLabel
Definition: cddefines.h:1235
t_struc::testr
realnum * testr
Definition: struc.h:25
t_tag_LineSv
Definition: lines.h:111
t_tag_LineSv::SumLine
double SumLine[4]
Definition: lines.h:125
t_tag_LineSv::chSumTyp
char chSumTyp
Definition: lines.h:114
t_tag_LineSv::chALab
char chALab[5]
Definition: lines.h:117
t_tag_LineSv::wavelength
realnum wavelength
Definition: lines.h:131
t_thermal::heating
double heating[LIMELM][LIMELM]
Definition: thermal.h:158
t_thermal::ctot
double ctot
Definition: thermal.h:112
t_thermal::dCooldT
double dCooldT
Definition: thermal.h:119
t_thermal::htot
double htot
Definition: thermal.h:149
t_thermal::totcol
double totcol
Definition: thermal.h:110
t_timesc::time_H2_Dest_here
double time_H2_Dest_here
Definition: timesc.h:37
t_timesc::sound_speed_adiabatic
double sound_speed_adiabatic
Definition: timesc.h:45
t_warnings::lgWarngs
bool lgWarngs
Definition: warnings.h:56
nSpecies
long int nSpecies
Definition: taulines.cpp:21
dBaseStates
vector< qList > dBaseStates
Definition: taulines.cpp:15
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
UTALines
TransitionList UTALines("UTALines", &AnonStates)
TauLine2
TransitionList TauLine2("TauLine2", &AnonStates)
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:17
HFLines
TransitionList HFLines("HFLines", &AnonStates)
nUTA
long int nUTA
Definition: taulines.cpp:26
nLevel1
long int nLevel1
Definition: taulines.cpp:28
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
TauLines
TransitionList TauLines("TauLines", &AnonStates)
dBaseSpecies
species * dBaseSpecies
Definition: taulines.cpp:14
TempChange
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:51
thermal
t_thermal thermal
Definition: thermal.cpp:5
timesc
t_timesc timesc
Definition: timesc.cpp:5
OccupationNumberLine
double OccupationNumberLine(const TransitionProxy &t)
Definition: transition.cpp:142
TexcLine
double TexcLine(const TransitionProxy &t)
Definition: transition.cpp:169
chLineLbl
char * chLineLbl(const TransitionProxy &t)
Definition: transition.cpp:237
warnings
t_warnings warnings
Definition: warnings.cpp:11
wind
Wind wind
Definition: wind.cpp:5
Generated by doxygen 1.9.5