cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
punch_do.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*PunchDo produce punch output during calculation,
4  * chTime is 'MIDL' during calculation, 'LAST' at the end */
5 /*PunchNewContinuum produce the 'punch new continuum' output */
6 /*PunchLineStuff punch optical depths or source functions for all transferred lines */
7 /*pun1Line called by PunchLineStuff to produce output for one line */
8 /*PunchNewContinuum produce the 'punch new continuum' output */
9 /*PunLineIntensity produce the 'punch lines intensity' output */
10 /* punch h emission, for AGN3 chapter 4, routine is below */
11 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */
12 /* the number of emission lines across one line of printout */
13 /*PunchSpecial generate output for the punch special command */
14 /*punResults punch results from punch results command */
15 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */
16 /*PunchGaunts called by punch gaunts command to output gaunt factors */
17 #include "cddefines.h"
18 #include "cddrive.h"
19 #include "physconst.h"
20 #include "mean.h"
21 #include "taulines.h"
22 #include "struc.h"
23 #include "iso.h"
24 #include "mole.h"
25 #include "hyperfine.h"
26 #include "rt.h"
27 #include "lines_service.h"
28 #include "doppvel.h"
29 #include "dense.h"
30 #include "h2.h"
31 #include "magnetic.h"
32 #include "hydrogenic.h"
33 #include "secondaries.h"
34 #include "grainvar.h"
35 #include "lines.h"
36 #include "dynamics.h"
37 #include "colden.h"
38 #include "continuum.h"
39 #include "ionbal.h"
40 #include "yield.h"
41 #include "prt.h"
42 #include "iterations.h"
43 #include "heavy.h"
44 #include "conv.h"
45 #include "geometry.h"
46 #include "called.h"
47 #include "helike.h"
48 #include "opacity.h"
49 #include "rfield.h"
50 #include "phycon.h"
51 #include "timesc.h"
52 #include "radius.h"
53 #include "atomfeii.h"
54 #include "assertresults.h"
55 #include "thermal.h"
56 #include "wind.h"
57 #include "hmi.h"
58 #include "pressure.h"
59 #include "elementnames.h"
60 #include "ipoint.h"
61 #include "gammas.h"
62 #include "atmdat.h"
63 #include "hcmap.h"
64 #include "input.h"
65 #include "punch.h"
66 #include "optimize.h"
67 #include "grid.h"
68 
69 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */
70 /* the number of emission lines across one line of printout */
72  FILE * ioPUN,
73  /* 4 char + null string */
74  const char *chLab,
75  realnum wl,
76  double xInten,
77  const char *chFunction);/* null terminated string saying what to do */
78 
79 /*PunchGaunts called by punch gaunts command to output gaunt factors */
80 STATIC void PunchGaunts(FILE* ioPUN);
81 
82 /*punResults punch results from punch results command */
83 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */
84 STATIC void punResults(FILE* ioPUN);
85 
87  FILE * ioPUN,
88  const char *chJob ,
89  realnum xLimit);
90 
91 /* punch h emission, for chapter 4, routine is below */
92 STATIC void AGN_Hemis(FILE *ioPUN );
93 
94 /*PunchNewContinuum produce the 'punch new continuum' output */
95 STATIC void PunchNewContinuum(FILE * ioPUN , long ipCon );
96 
97 /*PunLineIntensity produce the 'punch lines intensity' output */
98 STATIC void PunLineIntensity(FILE * ioPUN);
99 
100 char *chDummy;
101 
102 void PunchDo(
103  /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
104  const char *chTime)
105 {
106  bool lgOK,
107  lgTlkSav;
108  long int
109  ipPun,
110  i,
111  i1,
112  ion,
113  ipConMax,
114  ipHi,
115  ipLinMax,
116  ipLo,
117  ips,
118  j,
119  jj,
120  limit,
121  nd,
122  nelem;
123  double ConMax,
124  RateInter,
125  av,
126  conem,
127  eps,
128  flxatt,
129  flxcor,
130  flxin,
131  flxref,
132  flxtrn,
133  fout,
134  fref,
135  fsum,
136  opConSum,
137  opLinSum,
138  stage,
139  sum,
140  texc,
141  xLinMax;
142 
143  DEBUG_ENTRY( "PunchDo()" );
144 
145  /*
146  * the "last" option on punch command, to punch on last iteration,
147  * is parsed at the top of the loop in only one place.
148  * no further action is needed at all for punch last to work
149  * ok throughout this routine
150  */
151 
152  /*
153  * each branch can have a test whether chTime is or is not "LAST"
154  *
155  * if( strcmp(chTime,"LAST") == 0 ) <== print after iteration is complete
156  *
157  * if "LAST" then this is last call to routine after iteration complete
158  * punch only if "LAST" when results at end of iteration are needed
159  *
160  * if( strcmp(chTime,"LAST") != 0 ) <== print for every zone
161  *
162  * test for .not."LAST" is for every zone result, where you do not
163  * want to punch last zone twice
164  */
165 
166  /* return if no punch to do */
167  if( punch.npunch < 1 )
168  {
169  return;
170  }
171 
172  /* during a grid calculation this routine saves grid points after
173  * cloudy is called. we may output it below */
174  if( grid.lgGrid )
175  {
176  if( chTime[0]=='L' )
178 
179  /* save grid information */
180  GridGatherAfterCloudy( chTime );
181  }
182 
183  for( ipPun=0; ipPun < punch.npunch; ipPun++ )
184  {
185  /* this global variable to remember where in the punch stack we are */
186  punch.ipConPun = ipPun;
187 
188  /* iterations.lgLastIt is true if this is last iteration
189  * lgPunLstIter set true if 'last' key occurred on punch command
190  * normally is false. This will skip punching if last set and
191  * this is not last iteration */
192  if( iterations.lgLastIt || (!punch.lgPunLstIter[ipPun]) )
193  {
194 
195  if( strcmp(punch.chPunch[ipPun],"ABUN") == 0 )
196  {
197  /* punch abundances vs depth */
198  if( strcmp(chTime,"LAST") != 0 )
199  {
200  fprintf( punch.ipPnunit[ipPun], "%.2f",
202  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
203  {
204  /* >>chng 05 feb 03, protect against non-positive abundances,
205  * bug caught by Marcelo Castellanos */
206  fprintf( punch.ipPnunit[ipPun], "\t%.2f",
207  log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) );
208  }
209  fprintf( punch.ipPnunit[ipPun], "\n" );
210  }
211  }
212 
213  else if( strcmp(punch.chPunch[ipPun],"21CM") == 0 )
214  {
215  /* punch information about 21 cm line */
216  if( strcmp(chTime,"LAST") != 0 )
217  {
218  fprintf( punch.ipPnunit[ipPun],
219  "%.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",
220  /* depth, cm */
223  phycon.te ,
224  /* temperature from Lya - 21 cm optical depth ratio */
225  3.84e-7* Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon /
226  SDIV( HFLines[0].Emis->TauCon ),
227  /*TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ),*/
228  HFLines[0].Lo->Pop ,
229  HFLines[0].Hi->Pop ,
231  HFLines[0].Emis->TauCon ,
232  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon,
233  HFLines[0].Emis->PopOpc,
234  /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */
236  /* why was above multiplied by this following term? */
237  /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/
238  HFLines[0].Emis->TauIn,
241  /*>>chng 27 mar, GS, integrated 21cm spin temperature*/
244  -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower)
245  );
246  }
247  }
248 
249  else if( strcmp(punch.chPunch[ipPun],"AGES") == 0 )
250  {
251  /* punch timescales vs depth */
252  if( strcmp(chTime,"LAST") != 0 )
253  {
254  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
255  /* depth, cm */
257  /* cooling timescale */
259  /* H2 destruction timescale */
261  /* CO destruction timescale */
262  timesc.AgeCOMoleDest[findspecies("CO")->index],
263  /* OH destruction timescale */
264  timesc.AgeCOMoleDest[findspecies("OH")->index],
265  /* H recombination timescale */
266  1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) );
267  }
268  }
269 
270  else if( strcmp(punch.chPunch[ipPun]," AGN") == 0 )
271  {
272  if( strcmp(chTime,"LAST") == 0 )
273  {
274  if( strcmp( punch.chPunchArgs[ipPun], "HECS" ) == 0 )
275  {
276  /* this routine is in helike.c */
277  AGN_He1_CS(punch.ipPnunit[ipPun]);
278  }
279  if( strcmp( punch.chPunchArgs[ipPun], "HEMI" ) == 0 )
280  {
281  /* punch h emiss, for chapter 4, routine is below */
282  AGN_Hemis(punch.ipPnunit[ipPun]);
283  }
284  else
285  {
286  fprintf( ioQQQ, " PunchDo does not recognize flag %4.4s set for AGN punch. This is impossible.\n",
287  punch.chPunch[ipPun] );
288  ShowMe();
289  cdEXIT(EXIT_FAILURE);
290  }
291  }
292  }
293 
294  else if( strcmp(punch.chPunch[ipPun],"ASSE") == 0 )
295  {
296  if( strcmp(chTime,"LAST") == 0 )
297  {
298  /* punch the assert output */
299  lgCheckAsserts(punch.ipPnunit[ipPun]);
300  }
301  }
302 
303  else if( strcmp(punch.chPunch[ipPun],"AVER") == 0 )
304  {
305  if( strcmp(chTime,"LAST") == 0 )
306  {
307  /* punch the averages output */
308  punch_average( ipPun , "PUNCH", chDummy );
309  }
310  }
311 
312  else if( strncmp(punch.chPunch[ipPun],"CHA",3) == 0 )
313  {
314  if( strcmp(chTime,"LAST") == 0 )
315  {
316  /* one of the charge transfer options, all in chargtran.c */
317  ChargTranPun( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] );
318  }
319  }
320 
321  else if( strcmp(punch.chPunch[ipPun],"COOL") == 0 )
322  {
323  /* punch cooling, routine in file of same name */
324  if( strcmp(chTime,"LAST") != 0 )
325  CoolPunch(punch.ipPnunit[ipPun]);
326  }
327 
328  else if( strncmp(punch.chPunch[ipPun],"DYN" , 3) == 0 )
329  {
330  /* punch dynamics xxx, information about dynamical solutions */
331  if( strcmp(chTime,"LAST") != 0 )
332  DynaPunch(punch.ipPnunit[ipPun] ,punch.chPunch[ipPun][3] );
333  }
334 
335  else if( strcmp(punch.chPunch[ipPun],"ENTH") == 0 )
336  {
337  if( strcmp(chTime,"LAST") != 0 )
338  fprintf( punch.ipPnunit[ipPun],
339  "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
345  0.5*POW2(wind.windv)*dense.xMassDensity , /* KE */
346  5./2.*pressure.PresGasCurr , /* thermal plus PdV work */
347  magnetic.EnthalpyDensity); /* magnetic terms */
348  }
349 
350  else if( strcmp(punch.chPunch[ipPun],"COLU") == 0 )
351  {
352  /* punch column densities */
353  if( strcmp(chTime,"LAST") == 0 )
354  {
355  PrtColumns(punch.ipPnunit[ipPun]);
356  }
357  }
358 
359  else if( strcmp(punch.chPunch[ipPun],"COLS") == 0 )
360  {
361  /* punch some column densities */
362  if( strcmp(chTime,"LAST") == 0 )
363  {
364  char chHeader[2];
365  punch_colden(chHeader , punch.ipPnunit[ipPun] , "PUNS" );
366  }
367  }
368 
369  else if( strcmp(punch.chPunch[ipPun],"COMP") == 0 )
370  {
371  /* Compton energy exchange coefficients */
372  if( strcmp(chTime,"LAST") != 0 )
373  {
374  for( jj=0; jj<rfield.nflux; jj = jj + punch.ncPunchSkip)
375  {
376  fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e\n",
377  rfield.anu[jj], rfield.comup[jj]/rfield.widflx[jj],
378  rfield.comdn[jj]/rfield.widflx[jj] );
379  }
380  }
381  }
382 
383  /* punch continuum commands */
384  else if( strcmp(punch.chPunch[ipPun],"CON ") == 0 )
385  {
386  /* this is the usual "punch continuum" case */
387  /* >>chng 06 apr 03, add every option to do every zone */
388  /* if punch every is set then nPunchEveryZone must be positive
389  * was init to -1 */
390  bool lgPrintThis =false;
391  if( punch.lgPunchEveryZone[ipPun] )
392  {
393  /* this branch, every option is on line so want to print every n zone */
394  if( strcmp(chTime,"LAST") != 0 )
395  {
396  /* not last zone - print first and intermediate cases */
397  if( nzone==1 )
398  {
399  lgPrintThis = true;
400  }
401  else if( nzone%punch.nPunchEveryZone[ipPun]==0 )
402  {
403  lgPrintThis = true;
404  }
405  }
406  else
407  {
408  /* this is last zone, print only if did not trip on above */
409  if( nzone!=1 && nzone%punch.nPunchEveryZone[ipPun]!=0 )
410  {
411  lgPrintThis = true;
412  }
413  }
414  }
415  else
416  {
417  /* this branch, not "every", so only print the last zone */
418  if( strcmp(chTime,"LAST") == 0 )
419  lgPrintThis = true;
420  }
421  ASSERT( !punch.lgPunchEveryZone[ipPun] || punch.nPunchEveryZone[ipPun]>0 );
422  if( lgPrintThis )
423  {
424  if( punch.lgPunchEveryZone[ipPun] && nzone!=1)
425  fprintf( punch.ipPnunit[ipPun], "%s\n",
427 
428  for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
429  {
430  /* four continua predicted here;
431  * incident, attenuated incident, emitted,
432  * then attenuated incident + emitted, last reflected
433  * reflected continuum is stored relative to inner radius
434  * others are stored for this radius */
435 
436  /* NB this code also used in punch emitted,
437  * transmitted continuum commands */
438 
439  /* the incident continuum */
440  flxin = rfield.flux_total_incident[j]*rfield.anu2[j]*
441  EN1RYD/rfield.widflx[j];
442 
443  /* the reflected continuum */
444  flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/rfield.widflx[j] +
446 
447  /* the attenuated incident continuum */
448  flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/
450 
451  /* the outward emitted continuum */
452  conem = ((rfield.ConEmitOut[j])/
456 
457  /* sum of emitted and transmitted continua */
458  flxtrn = conem + flxatt;
459 
460  /* photon energy in appropriate energy or wavelength units */
461  fprintf( punch.ipPnunit[ipPun],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) );
462  /* incident continuum */
463  fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxin );
464  /* trans cont */
465  fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxatt );
466  /* DiffOut cont */
467  fprintf( punch.ipPnunit[ipPun],"%.3e\t", conem );
468  /* net trans cont */
469  fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn );
470  /* reflected cont */
471  fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref );
472  /* total cont */
473  fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref + flxtrn );
474  fprintf( punch.ipPnunit[ipPun], "%4.4s\t%4.4s\t",
475  /* line label */
476  rfield.chLineLabel[j] ,
477  /* cont label*/
478  rfield.chContLabel[j] );
479  /* number of lines within that cell over cell width
480  * punch raw continuum has number by itself */
481  fprintf( punch.ipPnunit[ipPun], "%.2f\n", rfield.line_count[j]/rfield.widflx[j]*rfield.anu[j] );
482  }
483  }
484  }
485 
486  /* this is the punch spectrum command,
487  * the new continuum command that will replace the previous one */
488  else if( strcmp(punch.chPunch[ipPun],"CONN") == 0 )
489  {
490  if( strcmp(chTime,"LAST") == 0 )
491  /* io unit and which new continuum this is (was set when punch read in */
492  PunchNewContinuum( punch.ipPnunit[ipPun] , (long)punch.punarg[ipPun][0]);
493  }
494 
495  else if( strcmp(punch.chPunch[ipPun],"CONC") == 0 )
496  {
497  /* punch incident continuum */
498  /* set pointer for possible change in units of energy in continuum
499  * AnuUnit will give anu in whatever units were set with punch units */
500  if( strcmp(chTime,"LAST") == 0 )
501  {
502  /* incident continuum */
503  for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
504  {
505  flxin = rfield.flux_total_incident[j]*rfield.anu2[j]*
506  EN1RYD/rfield.widflx[j];
507  /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
508  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\n",
509  AnuUnit(rfield.AnuOrg[j]), flxin );
510  }
511  }
512  }
513 
514  else if( strcmp(punch.chPunch[ipPun],"CONG") == 0 )
515  {
516  /* punch emitted grain continuum in optically thin limit */
517  if( strcmp(chTime,"LAST") == 0 )
518  {
519  /* GrainMakeDiffuse broke out emission into types
520  * according to matType */
521  for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
522  {
523  fsum = gv.GraphiteEmission[j]*rfield.anu2[j]*
525 
526  fout = gv.SilicateEmission[j]*rfield.anu2[j]*
528 
529  /* anu is .5e format to resolve energy mesh near 1 Ryd
530  * AnuUnit gives anu in whatever units were set with units option */
531  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n",
532  AnuUnit(rfield.AnuOrg[j]) , fsum , fout ,
533  /* total emission per unit volume defined in GrainMakeDiffuse
534  * used in RT_diffuse to form total grain emission */
535  gv.GrainEmission[j]*rfield.anu2[j]*
537  }
538  }
539  }
540 
541  else if( strcmp(punch.chPunch[ipPun],"CONR") == 0 )
542  {
543  /* punch reflected continuum */
544  /* set pointer for possible change in units of energy in continuum
545  * AnuUnit will give anu in whatever units were set with punch units */
546  if( strcmp(chTime,"LAST") == 0 )
547  {
548  if( geometry.lgSphere )
549  {
550  fprintf( punch.ipPnunit[ipPun], " Reflected continuum not predicted when SPHERE is set.\n" );
551  fprintf( ioQQQ ,
552  "\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" );
553  cdEXIT(EXIT_FAILURE);
554  }
555 
556  for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
557  {
558  /* reflected continuum */
559  flxref = rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/
560  rfield.widflx[j]*EN1RYD;
561  /* reflected lines */
562  fref = rfield.anu[j]*punch.PunchLWidth*
563  rfield.reflin[j]*EN1RYD;
564  /* ratio of reflected to incident continuum, the albedo */
565  if( rfield.flux_total_incident[j] > 1e-25 )
566  {
568  }
569  else
570  {
571  av = 0.;
572  }
573  /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
574  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4s\n",
575  AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref,
576  av, rfield.chContLabel[j] );
577  }
578  }
579  }
580 
581  else if( strcmp(punch.chPunch[ipPun],"CNVE") == 0 )
582  {
583  /* the punch convergence error command */
584  if( strcmp(chTime,"LAST") != 0 )
585  {
586  fprintf( punch.ipPnunit[ipPun],
587  "%.5e\t%li\t%.4e\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n",
593  dense.EdenTrue,
594  dense.eden,
596  thermal.htot,
597  thermal.ctot,
598  (thermal.htot - thermal.ctot)*100./thermal.htot );
599  }
600  }
601 
602  else if( strcmp(punch.chPunch[ipPun],"CONB") == 0 )
603  {
604  /* punch continuum bins binning */
605  /* set pointer for possible change in units of energy in continuum
606  * AnuUnit will give anu in whatever units were set with punch units */
607  if( strcmp(chTime,"LAST") != 0 )
608  {
609  for( j=0; j<rfield.nupper; j = j + punch.ncPunchSkip)
610  {
611  fprintf( punch.ipPnunit[ipPun], "%14.5e %14.5e %14.5e\n",
612  AnuUnit(rfield.AnuOrg[j]), rfield.anu[j], rfield.widflx[j] );
613  }
614  }
615  }
616 
617  else if( strcmp(punch.chPunch[ipPun],"COND") == 0 )
618  {
619  /* punch diffuse continuum the local thermal emission */
620  /* set pointer for possible change in units of energy in continuum
621  * AnuUnit will give anu in whatever units were set with punch units */
622  if( strcmp(chTime,"LAST") == 0 )
623  {
624  /* this option to punch diffuse continuum */
625  for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
626  {
627  /* >>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */
628  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\n",
629  AnuUnit(rfield.AnuOrg[j]),
631  }
632  }
633  }
634 
635  else if( strcmp(punch.chPunch[ipPun],"CONE") == 0 )
636  {
637  /* punch emitted continuum */
638  /* set pointer for possible change in units of energy in continuum
639  * AnuUnit will give anu in whatever units were set with punch units */
640  if( strcmp(chTime,"LAST") == 0 )
641  {
642  /* punch emitted continuum */
643  for( j=0; j<rfield.nflux;j = j +punch.ncPunchSkip)
644  {
645  /* this is the reflected component */
646  flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/
648  rfield.reflin[j])*EN1RYD;
649 
650  /* this is the total emission in the outward direction */
651  conem = (rfield.ConEmitOut[j])/rfield.widflx[j]*rfield.anu2[j] +
653 
655 
656  /* output: photon energy, reflected, outward, total emission
657  * >>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */
658  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n",
659  AnuUnit(rfield.AnuOrg[j]),
660  flxref,
661  conem,
662  flxref + conem,
663  rfield.chLineLabel[j],
665  );
666  }
667  }
668  }
669 
670  /* punch fine continuum command */
671  else if( strcmp(punch.chPunch[ipPun],"CONf") == 0 )
672  {
673  if( strcmp(chTime,"LAST") == 0 )
674  {
675  long nu_hi , nskip;
676  if( punch.punarg[ipPun][0] > 0. )
677  /* possible lower bounds to energy range -
678  * 0 if not set with range option*/
679  j = ipFineCont( punch.punarg[ipPun][0] );
680  else
681  j = 0;
682 
683  /* upper limit set with range option */
684  if( punch.punarg[ipPun][1]> 0. )
685  nu_hi = ipFineCont( punch.punarg[ipPun][1]);
686  else
687  nu_hi = rfield.nfine;
688 
689  /* number of cells to bring together, default is 10 */
690  nskip = (long)punch.punarg[ipPun][2];
691 
692  do
693  {
694  realnum sum1 = rfield.fine_opt_depth[j];
695  realnum xnu = rfield.fine_anu[j];
696  for( jj=1; jj<nskip; ++jj )
697  {
698  xnu += rfield.fine_anu[j+jj];
699  sum1 += rfield.fine_opt_depth[j+jj];
700  }
701  fprintf( punch.ipPnunit[ipPun],
702  "%.6e\t%.3e\n",
703  AnuUnit(xnu/nskip),
704  sexp(sum1/nskip) );
705  j += nskip;
706  } while( j < nu_hi );
707  }
708  }
709 
710  else if( strcmp(punch.chPunch[ipPun],"CONi") == 0 )
711  {
712  /* punch continuum interactions */
713  /* set pointer for possible change in units of energy in continuum
714  * AnuUnit will give anu in whatever units were set with punch units */
715 
716  /* continuum interactions */
717  if( strcmp(chTime,"LAST") != 0 )
718  {
719  /* this is option to set lowest energy */
720  if( punch.punarg[ipPun][0] <= 0. )
721  {
722  i1 = 1;
723  }
724  else if( punch.punarg[ipPun][0] < 100. )
725  {
726  i1 = ipoint(punch.punarg[ipPun][0]);
727  }
728  else
729  {
730  i1 = (long int)punch.punarg[ipPun][0];
731  }
732 
733  fref = 0.;
734  fout = 0.;
735  fsum = 0.;
736  sum = 0.;
737  flxin = 0.;
738 
739  for( j=i1-1; j < rfield.nflux; j++ )
740  {
741  fref += rfield.flux[j]*opac.opacity_abs[j];
742  fout += rfield.otslin[j]*opac.opacity_abs[j];
743  fsum += rfield.otscon[j]*opac.opacity_abs[j];
744  sum += rfield.ConInterOut[j]*opac.opacity_abs[j];
745  flxin += (rfield.outlin[j] + rfield.outlin_noplot[j])*opac.opacity_abs[j];
746  }
747  fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e%10.2e%10.2e\n",
748  fref, fout, fsum, sum, flxin );
749  }
750  }
751 
752  else if( strcmp(punch.chPunch[ipPun],"CONI") == 0 )
753  {
754  /* punch ionizing continuum */
755  /* set pointer for possible change in units of energy in continuum
756  * AnuUnit will give anu in whatever units were set with punch units */
757 
758  if( (punch.punarg[ipPun][2]>0) || (strcmp(chTime,"LAST") == 0) )
759  {
760  /* this flag will remember whether we have ever printed anything */
761  bool lgPrt=false;
762  if( punch.punarg[ipPun][2]>0 )
763  fprintf(punch.ipPnunit[ipPun],"#punch every zone %li\n", nzone);
764 
765  /* punch ionizing continuum command
766  * this is option to set lowest energy,
767  * if no number was entered then this was zero */
768  if( punch.punarg[ipPun][0] <= 0. )
769  {
770  i1 = 1;
771  }
772  else if( punch.punarg[ipPun][0] < 100. )
773  {
774  i1 = ipoint(punch.punarg[ipPun][0]);
775  }
776  else
777  {
778  i1 = (long int)punch.punarg[ipPun][0];
779  }
780 
781  sum = 0.;
782  for( j=i1-1; j < rfield.nflux; j++ )
783  {
784  flxcor = rfield.flux[j] +
785  rfield.otslin[j] +
786  rfield.otscon[j] +
787  rfield.ConInterOut[j] +
789 
790  sum += flxcor*opac.opacity_abs[j];
791  }
792 
793  if( sum > 0. )
794  {
795  sum = 1./sum;
796  }
797  else
798  {
799  sum = 1.;
800  }
801 
802  fsum = 0.;
803 
804  for( j=i1-1; j<rfield.nflux; ++j)
805  {
806  flxcor = rfield.flux[j] +
807  rfield.otslin[j] +
808  rfield.otscon[j] +
809  rfield.ConInterOut[j]+
811 
812  fsum += flxcor*opac.opacity_abs[j];
813 
814  /* punched quantities are freq, flux, flux*cross sec,
815  * fraction of total, integral fraction of total */
816  RateInter = flxcor*opac.opacity_abs[j]*sum;
817 
818  /* punage(ipPun,2) is lowest interaction rate to consider, def=0.01 (1 percent) */
819  /* >>chng 01 nov 22, format to c-friendly */
820  if( (RateInter >= punch.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) )
821  {
822  lgPrt = true;
823  /* >>chng 96 oct 22, format of anu to 11.5 to resolve energy mesh near 1 Ryd */
824  fprintf( punch.ipPnunit[ipPun],
825  "%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",
826  j,
827  AnuUnit(rfield.AnuOrg[j]),
828  flxcor,
829  flxcor*opac.opacity_abs[j],
830  rfield.flux[j]/flxcor,
831  rfield.otslin[j]/flxcor,
832  rfield.otscon[j]/flxcor,
833  (rfield.outlin[j] + rfield.outlin_noplot[j])/flxcor,
834  rfield.ConInterOut[j]/flxcor,
835  RateInter,
836  fsum*sum,
837  rfield.chLineLabel[j],
838  rfield.chContLabel[j] );
839  }
840  }
841  if( !lgPrt )
842  {
843  /* entered logical block but did not print anything */
844  fprintf(punch.ipPnunit[ipPun],
845  " punchdo, the PUNCH IONIZING CONTINUUM command "
846  "did not find a strong point, sum and fsum were %.2e %.2e\n",
847  sum,fsum);
848  fprintf(punch.ipPnunit[ipPun],
849  " punchdo, the low-frequency energy was %.5e Ryd\n",
850  rfield.anu[i1-1]);
851  fprintf(punch.ipPnunit[ipPun],
852  " You can reset the threshold for the lowest fractional "
853  "interaction to print with the second number of the punch command\n"
854  " The fraction was %.3f and this was too large.\n",
855  punch.punarg[ipPun][1]);
856  }
857  }
858  }
859 
860  else if( strcmp(punch.chPunch[ipPun],"CORA") == 0 )
861  {
862  /* punch raw continuum */
863  /* set pointer for possible change in units of energy in continuum
864  * AnuUnit will give anu in whatever units were set with punch units */
865 
866  if( strcmp(chTime,"LAST") == 0 )
867  {
868  /* this option to punch all raw ionizing continuum */
869  for( j=0;j<rfield.nflux;j = j + punch.ncPunchSkip)
870  {
871  fprintf( punch.ipPnunit[ipPun],
872  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t",
873  AnuUnit(rfield.AnuOrg[j]),
874  rfield.flux[j],
875  rfield.otslin[j],
876  rfield.otscon[j],
877  rfield.ConRefIncid[j],
878  rfield.ConEmitReflec[j],
879  rfield.ConInterOut[j],
881  rfield.ConEmitOut[j] ,
882  rfield.chLineLabel[j],
884  );
885  /* number of lines within that cell */
886  fprintf( punch.ipPnunit[ipPun], "%li\n", rfield.line_count[j] );
887  }
888  }
889  }
890 
891  else if( strcmp(punch.chPunch[ipPun],"CONo") == 0 )
892  {
893  /* punch outward local continuum */
894  /* set pointer for possible change in units of energy in continuum
895  * AnuUnit will give anu in whatever units were set with punch units */
896 
897  if( strcmp(chTime,"LAST") == 0 )
898  {
899  /* option to punch out outward continuum here */
900  for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip)
901  {
902  fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e\n",
903  AnuUnit(rfield.AnuOrg[j]),
905  (rfield.flux[j] + rfield.otscon[j] + rfield.otslin[j] +
907  rfield.anu[j] );
908  }
909  }
910  }
911 
912  else if( strcmp(punch.chPunch[ipPun],"CONO") == 0 )
913  {
914  /* punch outward continuum */
915  /* set pointer for possible change in units of energy in continuum
916  * AnuUnit will give anu in whatever units were set with punch units */
917 
918  if( strcmp(chTime,"LAST") == 0 )
919  {
920  /* option to punch out outward continuum */
921  for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
922  {
923  fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e%10.2e%10.2e\n",
924  AnuUnit(rfield.AnuOrg[j]),
930  rfield.anu[j])*radius.r1r0sq*EN1RYD );
931  }
932  }
933  }
934 
935  else if( strcmp(punch.chPunch[ipPun],"CONT") == 0 )
936  {
937  /* punch transmitted continuum - this is not the main "punch continuum"
938  * command - search on "CON " above
939  * set pointer for possible change in units of energy in continuum
940  * AnuUnit will give anu in whatever units were set with punch units */
941 
942  if( strcmp(chTime,"LAST") == 0 )
943  {
944  /* this option to punch transmitted continuum */
945  for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip)
946  {
947  /* attenuated incident continuum
948  * >>chng 97 jul 10, remove PunchLWidth from this one only since
949  * we must conserve energy even in lines
950  * >>chng 07 apr 26 include transmission coefficient */
951  flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/
953 
954  /*conem = (rfield.ConOutNoInter[j] + rfield.ConInterOut[j]+rfield.outlin[j])*
955  rfield.anu2[j];
956  conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;*/
957  /* >>chng 00 jan 03, above did not include all contributors.
958  * Pasted in below from usual
959  * punch continuum command */
960  /* >>chng 04 jul 15, removed factor of punch.PunchLWidth -
961  * this should not be there to conserve energy, as explained in hazy
962  * where command was documented, and in comment above. caught by PvH */
963  /* >>chng 04 jul 23, incorrect use of outlin - before multiplied by an2,
964  * quantity should be photons per Ryd, since init quantity is
965  * photons per cell. Must div by widflx. caught by PvH */
966  /*conem = (rfield.ConEmitOut[j]/rfield.widflx[j]*rfield.anu2[j] +
967  rfield.outlin[j]*rfield.anu[j])*radius.r1r0sq*
968  EN1RYD*geometry.covgeo;*/
969  conem = (rfield.ConEmitOut[j] + rfield.outlin[j]) / rfield.widflx[j]*
971 
972  flxtrn = conem + flxatt;
973 
974  /* use AnuOrg here instead of anu since probably
975  * going to be used by table read
976  * and we want good anu array for sanity check*/
977  fprintf( punch.ipPnunit[ipPun],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) );
978  fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn);
979  fprintf( punch.ipPnunit[ipPun],"%.3e\n", rfield.trans_coef_total[j] );
980  }
981  }
982  }
983 
984  else if( strcmp(punch.chPunch[ipPun],"CON2") == 0 )
985  {
986  /* punch total two-pohoton continuum */
987  if( strcmp(chTime,"LAST") == 0 )
988  {
989  /* this option to punch diffuse continuum */
990  for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
991  {
992  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\n",
993  AnuUnit(rfield.AnuOrg[j]),
996  }
997  }
998  }
999 
1000  else if( strcmp(punch.chPunch[ipPun],"DUSE") == 0 )
1001  {
1002  /* punch grain extinction - includes only grain opacity, not total */
1003  if( strcmp(chTime,"LAST") != 0 )
1004  {
1005  fprintf( punch.ipPnunit[ipPun], " %.5e\t",
1007 
1008  /* visual extinction of an extended source (like a PDR)*/
1009  fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
1010 
1011  /* visual extinction of point source (star)*/
1012  fprintf( punch.ipPnunit[ipPun], "%.2e\n" , rfield.extin_mag_V_point);
1013  }
1014  }
1015 
1016  else if( strcmp(punch.chPunch[ipPun],"DUSO") == 0 )
1017  {
1018  /* punch grain opacity */
1019  if( strcmp(chTime,"LAST") == 0 )
1020  {
1021  for( j=0; j < rfield.nflux; j++ )
1022  {
1023  double scat;
1024  fprintf( punch.ipPnunit[ipPun],
1025  "%.5e\t%.2e\t%.2e\t%.2e\t",
1026  /* photon energy or wavelength */
1027  AnuUnit(rfield.AnuOrg[j]),
1028  /* total opacity, discount forward scattering */
1029  gv.dstab[j] + gv.dstsc[j],
1030  /* absorption opacity */
1031  gv.dstab[j],
1032  /* scatter, with forward discounted */
1033  gv.dstsc[j] );
1034  /* add together total scattering, discounting 1-g */
1035  scat = 0.;
1036  /* sum over all grain species */
1037  for( nd=0; nd < gv.nBin; nd++ )
1038  {
1039  scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund;
1040  }
1041  /* finally, scattering including effects of forward scattering */
1042  fprintf( punch.ipPnunit[ipPun],
1043  "%.2e", scat );
1044  fprintf( punch.ipPnunit[ipPun],
1045  "%.2e\n", gv.dstsc[j]/(gv.dstab[j] + gv.dstsc[j]) );
1046  }
1047  }
1048  }
1049 
1050  /* punch grain abundance and punch grain D/G ratio commands */
1051  else if( strcmp(punch.chPunch[ipPun],"DUSA") == 0 ||
1052  strcmp(punch.chPunch[ipPun],"DUSD") == 0 )
1053  {
1054  bool lgDGRatio = ( strcmp(punch.chPunch[ipPun],"DUSD") == 0 );
1055 
1056  /* grain abundance */
1057  if( strcmp(chTime,"LAST") != 0 )
1058  {
1059  /* used to print header exactly one time */
1060  static bool lgMustPrtHeaderDRRatio = true,
1061  lgMustPrtHeaderAbundance=true;
1062  /* print grain headder first if this has not yet been done */
1063  if( ( lgMustPrtHeaderDRRatio && lgDGRatio ) ||
1064  ( lgMustPrtHeaderAbundance && !lgDGRatio ) )
1065  {
1066  /* only print one header for each case, but must
1067  * track separately if both used in same sim */
1068  if( lgMustPrtHeaderDRRatio && lgDGRatio )
1069  lgMustPrtHeaderDRRatio = false;
1070  else if( lgMustPrtHeaderAbundance &&!lgDGRatio )
1071  lgMustPrtHeaderAbundance = false;
1072 
1073  fprintf( punch.ipPnunit[ipPun], "#Depth" );
1074  for( nd=0; nd < gv.nBin; ++nd )
1075  fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1076  fprintf( punch.ipPnunit[ipPun], "\ttotal\n" );
1077 
1078  /* now print grain radius converting from cm to microns */
1079  fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
1080  for( nd=0; nd < gv.nBin; ++nd )
1081  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1082  fprintf( punch.ipPnunit[ipPun], "\txxxx\n" );
1083  }
1084  fprintf( punch.ipPnunit[ipPun], " %.5e",
1086  /* grain abundance per bin in g/cm^3 */
1087  double total = 0.;
1088  for( nd=0; nd < gv.nBin; ++nd )
1089  {
1090  double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]*
1091  gv.bin[nd]->cnv_H_pCM3;
1092  if( lgDGRatio )
1093  abund /= dense.xMassDensity;
1094  fprintf( punch.ipPnunit[ipPun], "\t%.3e", abund );
1095  total += abund;
1096  }
1097  fprintf( punch.ipPnunit[ipPun], "\t%.3e\n", total );
1098  }
1099  }
1100 
1101  else if( strcmp(punch.chPunch[ipPun],"DUSP") == 0 )
1102  {
1103  /* grain potential */
1104  if( strcmp(chTime,"LAST") != 0 )
1105  {
1106  /* used to print header exactly one time */
1107  static bool lgMustPrtHeader = true;
1108  /* do labels first if this is first zone */
1109  if( lgMustPrtHeader )
1110  {
1111  /* first print string giving grain id */
1112  fprintf( punch.ipPnunit[ipPun], "#Depth" );
1113  for( nd=0; nd < gv.nBin; ++nd )
1114  fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1115  fprintf( punch.ipPnunit[ipPun], "\n" );
1116 
1117  /* now print grain radius converting from cm to microns */
1118  fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
1119  for( nd=0; nd < gv.nBin; ++nd )
1120  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1121  fprintf( punch.ipPnunit[ipPun], "\n" );
1122 
1123  /* don't need to do this, ever again */
1124  lgMustPrtHeader = false;
1125  }
1126  fprintf( punch.ipPnunit[ipPun], " %.5e",
1128  /* grain potential in eV */
1129  for( nd=0; nd < gv.nBin; ++nd )
1130  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD );
1131  fprintf( punch.ipPnunit[ipPun], "\n" );
1132  }
1133  }
1134 
1135  else if( strcmp(punch.chPunch[ipPun],"DUSR") == 0 )
1136  {
1137  /* grain H2 formation rates */
1138  if( strcmp(chTime,"LAST") != 0 )
1139  {
1140  /* used to print header exactly one time */
1141  static bool lgMustPrtHeader = true;
1142  /* do labels first if this is first zone */
1143  if( lgMustPrtHeader )
1144  {
1145  /* first print string giving grain id */
1146  fprintf( punch.ipPnunit[ipPun], "#Depth" );
1147  for( nd=0; nd < gv.nBin; ++nd )
1148  fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1149  fprintf( punch.ipPnunit[ipPun], "\n" );
1150 
1151  /* now print grain radius converting from cm to microns */
1152  fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
1153  for( nd=0; nd < gv.nBin; ++nd )
1154  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1155  fprintf( punch.ipPnunit[ipPun], "\n" );
1156 
1157  /* don't need to do this, ever again */
1158  lgMustPrtHeader = false;
1159  }
1160  fprintf( punch.ipPnunit[ipPun], " %.5e",
1162  /* grain formation rate for H2 */
1163  for( nd=0; nd < gv.nBin; ++nd )
1164  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used );
1165  fprintf( punch.ipPnunit[ipPun], "\n" );
1166  }
1167  }
1168 
1169  else if( strcmp(punch.chPunch[ipPun],"DUST") == 0 )
1170  {
1171  /* grain temperatures - K*/
1172  if( strcmp(chTime,"LAST") != 0 )
1173  {
1174  /* used to print header exactly one time */
1175  static bool lgMustPrtHeader = true;
1176  /* do labels first if this is first zone */
1177  if( lgMustPrtHeader )
1178  {
1179  /* first print string giving grain id */
1180  fprintf( punch.ipPnunit[ipPun], "#Depth" );
1181  for( nd=0; nd < gv.nBin; ++nd )
1182  fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1183  fprintf( punch.ipPnunit[ipPun], "\n" );
1184 
1185  /* now print grain radius converting from cm to microns */
1186  fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
1187  for( nd=0; nd < gv.nBin; ++nd )
1188  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1189  fprintf( punch.ipPnunit[ipPun], "\n" );
1190 
1191  /* don't need to do this, ever again */
1192  lgMustPrtHeader = false;
1193  }
1194  fprintf( punch.ipPnunit[ipPun], " %.5e",
1196  for( nd=0; nd < gv.nBin; ++nd )
1197  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust );
1198  fprintf( punch.ipPnunit[ipPun], "\n" );
1199  }
1200  }
1201 
1202  else if( strcmp(punch.chPunch[ipPun],"DUSC") == 0 )
1203  {
1204  /* punch grain charge - eden from grains and
1205  * charge per grain in electrons / grain */
1206  if( strcmp(chTime,"LAST") != 0 )
1207  {
1208  /* used to print header exactly one time */
1209  static bool lgMustPrtHeader = true;
1210  /* do labels first if this is first zone */
1211  if( lgMustPrtHeader )
1212  {
1213  /* first print string giving grain id */
1214  fprintf( punch.ipPnunit[ipPun], "#Depth\tne(grn)" );
1215  for( nd=0; nd < gv.nBin; ++nd )
1216  fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1217  fprintf( punch.ipPnunit[ipPun], "\n" );
1218 
1219  /* now print grain radius converting from cm to microns */
1220  fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)\tne\t" );
1221  for( nd=0; nd < gv.nBin; ++nd )
1222  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1223  fprintf( punch.ipPnunit[ipPun], "\n" );
1224 
1225  /* don't need to do this, ever again */
1226  lgMustPrtHeader = false;
1227  }
1228 
1229  fprintf( punch.ipPnunit[ipPun], " %.5e\t%.4e",
1231  /* electron density contributed by grains, in e/cm^3,
1232  * positive number means grain supplied free electrons */
1233  gv.TotalEden );
1234 
1235  /* average charge per grain in electrons */
1236  for( nd=0; nd < gv.nBin; ++nd )
1237  {
1238  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ );
1239  }
1240  fprintf( punch.ipPnunit[ipPun], "\n" );
1241  }
1242  }
1243 
1244  else if( strcmp(punch.chPunch[ipPun],"DUSH") == 0 )
1245  {
1246  /* grain heating */
1247  if( strcmp(chTime,"LAST") != 0 )
1248  {
1249  /* used to print header exactly one time */
1250  static bool lgMustPrtHeader = true;
1251  /* punch grain charge, but do labels first if this is first zone */
1252  if( lgMustPrtHeader )
1253  {
1254  /* first print string giving grain id */
1255  fprintf( punch.ipPnunit[ipPun], "#Depth" );
1256  for( nd=0; nd < gv.nBin; ++nd )
1257  fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1258  fprintf( punch.ipPnunit[ipPun], "\n" );
1259 
1260  /* now print grain radius converting from cm to microns */
1261  fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
1262  for( nd=0; nd < gv.nBin; ++nd )
1263  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1264  fprintf( punch.ipPnunit[ipPun], "\n" );
1265 
1266  /* don't need to do this, ever again */
1267  lgMustPrtHeader = false;
1268  }
1269  fprintf( punch.ipPnunit[ipPun], " %.5e",
1271  /* grain heating */
1272  for( nd=0; nd < gv.nBin; ++nd )
1273  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl );
1274  fprintf( punch.ipPnunit[ipPun], "\n" );
1275  }
1276  }
1277 
1278  else if( strcmp(punch.chPunch[ipPun],"DUSV") == 0 )
1279  {
1280  /* grain drift velocities */
1281  if( strcmp(chTime,"LAST") != 0 )
1282  {
1283  /* used to print header exactly one time */
1284  static bool lgMustPrtHeader = true;
1285  /* punch grain velocity, but do labels first if this is first zone */
1286  if( lgMustPrtHeader )
1287  {
1288  /* first print string giving grain id */
1289  fprintf( punch.ipPnunit[ipPun], "#Depth" );
1290  for( nd=0; nd < gv.nBin; ++nd )
1291  fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1292  fprintf( punch.ipPnunit[ipPun], "\n" );
1293 
1294  /* now print grain radius converting from cm to microns */
1295  fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
1296  for( nd=0; nd < gv.nBin; ++nd )
1297  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1298  fprintf( punch.ipPnunit[ipPun], "\n" );
1299 
1300  /* don't need to do this, ever again */
1301  lgMustPrtHeader = false;
1302  }
1303  fprintf( punch.ipPnunit[ipPun], " %.5e",
1305  /* grain drift velocity in km/s */
1306  for( nd=0; nd < gv.nBin; ++nd )
1307  fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 );
1308  fprintf( punch.ipPnunit[ipPun], "\n" );
1309  }
1310  }
1311 
1312  /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
1313  else if( strcmp(punch.chPunch[ipPun],"DUSQ") == 0 )
1314  {
1315  /* punch grain Qs */
1316  if( strcmp(chTime,"LAST") == 0 )
1317  {
1318  for( j=0; j < rfield.nflux; j++ )
1319  {
1320  fprintf( punch.ipPnunit[ipPun], "%11.4e",
1321  rfield.anu[j] );
1322  for( nd=0; nd < gv.nBin; nd++ )
1323  {
1324  fprintf( punch.ipPnunit[ipPun], "%9.1e%9.1e",
1325  gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea,
1326  gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea );
1327  }
1328  fprintf( punch.ipPnunit[ipPun], "\n" );
1329  }
1330  }
1331  }
1332 
1333  else if( strcmp(punch.chPunch[ipPun],"ELEM") == 0 )
1334  {
1335  if( strcmp(chTime,"LAST") != 0 )
1336  {
1337  realnum renorm = 1.f;
1338 
1339  /* this is the index for the atomic number on the physical scale */
1340  /* >>chng 04 nov 23, use c scale throughout */
1341  nelem = (long int)punch.punarg[ipPun][0];
1342  ASSERT( nelem >= ipHYDROGEN );
1343 
1344  /* don't do this if element is not turned on */
1345  if( dense.lgElmtOn[nelem] )
1346  {
1347  /* >>chng 04 nov 23, add density option, leave as cm-3
1348  * default is still norm to total of that element */
1349  if( punch.punarg[ipPun][1] == 0 )
1350  renorm = dense.gas_phase[nelem];
1351 
1352  fprintf( punch.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone );
1353 
1354  for( j=0; j <= (nelem + 1); ++j)
1355  {
1356  fprintf( punch.ipPnunit[ipPun], "\t%.2e",
1357  dense.xIonDense[nelem][j]/renorm );
1358  }
1359  if( nelem==ipHYDROGEN )
1360  {
1361  /* H2 */
1362  fprintf( punch.ipPnunit[ipPun], "\t%.2e",
1363  hmi.H2_total/renorm );
1364  }
1365  /* >>chng 04 nov 23 add C and O fine structure pops */
1366  else if( nelem==ipCARBON )
1367  {
1368  fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
1369  colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm);
1370  fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e",
1371  colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm);
1372  fprintf( punch.ipPnunit[ipPun], "\t%.2e",
1373  findspecies("CO")->hevmol/renorm );
1374  }
1375  else if( nelem==ipOXYGEN )
1376  {
1377  fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
1378  colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm);
1379  }
1380  fprintf( punch.ipPnunit[ipPun], "\n" );
1381  }
1382  }
1383  }
1384 
1385  else if( strcmp(punch.chPunch[ipPun],"RECA") == 0 )
1386  {
1387  /* this will create table for AGN3 then exit,
1388  * routine is in makerecom.c */
1389  ion_recombAGN( punch.ipPnunit[ipPun] );
1390  cdEXIT(EXIT_FAILURE);
1391  }
1392 
1393  else if( strcmp(punch.chPunch[ipPun],"RECE") == 0 )
1394  {
1395  /* punch recombination efficiencies,
1396  * option turned on with the "punch recombination efficiencies" command
1397  * output for the punch recombination coefficients command is actually
1398  * produced by a series of routines, as they generate the recombination
1399  * coefficients. these include
1400  * dielsupres, helium, hydrorecom, iibod, and makerecomb*/
1401  fprintf( punch.ipPnunit[ipPun],
1402  "%12.4e %12.4e %12.4e %12.4e\n",
1405  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecRad],
1406  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecNetEsc]);
1407  }
1408 
1409  else if( strcmp(punch.chPunch[ipPun],"FEco") == 0 )
1410  {
1411  /* FeII continuum, check that FeII is turned on */
1412  if( strcmp(chTime,"LAST") == 0 )
1413  {
1414  for( j=0; j < nFeIIConBins; j++ )
1415  {
1416  /* emission from large FeII atom, integrated over continuum intervals
1417  * units are always in units, erg cm-2 s-1 as in intensity case,
1418  * even if luminosity case is considered
1419  * [j][0] is flux, [j][1] and [j][2] are upper and
1420  * lower bounds in Angstroms.
1421  * these are set in FeIIZero */
1424  fprintf( punch.ipPnunit[ipPun], "%.2f\t%e \n",
1425  (FeII_Cont[j][1]+FeII_Cont[j][2])/2.,
1426  FeII_Cont[j][0] );
1427  }
1428  }
1429  }
1430 
1431  /* punch column densities */
1432  else if( strcmp(punch.chPunch[ipPun],"FEcl") == 0 )
1433  {
1434  if( strcmp(chTime,"LAST") == 0 )
1435  {
1436  /* punch FeII level energies and stat weights, followed by column density */
1437  FeIIPunchColden( punch.ipPnunit[ipPun] );
1438  }
1439  }
1440 
1441  else if( strcmp(punch.chPunch[ipPun],"FE2l") == 0 )
1442  {
1443  if( strcmp(chTime,"LAST") == 0 )
1444  {
1445  /* punch FeII level energies and stat weights */
1446  FeIIPunchLevels( punch.ipPnunit[ipPun] );
1447  }
1448  }
1449 
1450  else if( strcmp(punch.chPunch[ipPun],"FEli") == 0 )
1451  {
1452  if( strcmp(chTime,"LAST") == 0 )
1453  {
1454  /* punch line intensities, routine is in atom_feii.c */
1455  FeIIPunchLines( punch.ipPnunit[ipPun] );
1456  }
1457  }
1458 
1459  else if( strcmp(punch.chPunch[ipPun],"FE2o") == 0 )
1460  {
1461  if( strcmp(chTime,"LAST") == 0 )
1462  {
1463  /* punch line optical depths, routine is in atom_feii.c */
1465  }
1466  }
1467 
1468  else if( strcmp(punch.chPunch[ipPun],"FRED") == 0 )
1469  {
1470  /* set with punch Fred command, this punches some stuff from
1471  * Fred Hamann's dynamics project */
1472  if( strcmp(chTime,"LAST") != 0 )
1473  {
1474  /* Fred's list */
1475  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\t%.3e\t%.3e\t%.3e"
1476  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
1477  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
1478  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
1479  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
1483  wind.fmul ,
1484  mean.xIonMeans[0][ipHYDROGEN][0] , mean.xIonMeans[0][ipHYDROGEN][1] ,
1485  mean.xIonMeans[0][ipHELIUM][0] , mean.xIonMeans[0][ipHELIUM][1] ,
1486  mean.xIonMeans[0][ipHELIUM][2] ,
1487  mean.xIonMeans[0][ipCARBON][1] , mean.xIonMeans[0][ipCARBON][2] ,
1488  mean.xIonMeans[0][ipCARBON][3] ,
1489  mean.xIonMeans[0][ipOXYGEN][0] , mean.xIonMeans[0][ipOXYGEN][1] ,
1490  mean.xIonMeans[0][ipOXYGEN][2] , mean.xIonMeans[0][ipOXYGEN][3] ,
1491  mean.xIonMeans[0][ipOXYGEN][4] , mean.xIonMeans[0][ipOXYGEN][5] ,
1492  mean.xIonMeans[0][ipOXYGEN][6] , mean.xIonMeans[0][ipOXYGEN][7] ,
1495  dense.xIonDense[ipHELIUM][2] ,
1497  dense.xIonDense[ipCARBON][3] ,
1503  }
1504  }
1505 
1506  else if( strcmp(punch.chPunch[ipPun],"FE2d") == 0 )
1507  {
1508  /* punch some departure coefficients for large FeII atom */
1509  if( strcmp(chTime,"LAST") != 0 )
1510  FeIIPunDepart(punch.ipPnunit[ipPun],false);
1511  }
1512 
1513  else if( strcmp(punch.chPunch[ipPun],"FE2D") == 0 )
1514  {
1515  /* punch all departure coefficients for large FeII atom */
1516  if( strcmp(chTime,"LAST") != 0 )
1517  FeIIPunDepart(punch.ipPnunit[ipPun],true);
1518  }
1519 
1520  else if( strcmp(punch.chPunch[ipPun],"FE2p") == 0 )
1521  {
1522  bool lgFlag = false;
1523  if( punch.punarg[ipPun][2] )
1524  lgFlag = true;
1525  /* punch small subset of level populations for large FeII atom */
1526  if( strcmp(chTime,"LAST") != 0 )
1527  FeIIPunPop(punch.ipPnunit[ipPun],false,0,0,
1528  lgFlag);
1529  }
1530 
1531  else if( strcmp(punch.chPunch[ipPun],"FE2P") == 0 )
1532  {
1533  bool lgFlag = false;
1534  if( punch.punarg[ipPun][2] )
1535  lgFlag = true;
1536  /* punch range of level populations for large FeII atom */
1537  if( strcmp(chTime,"LAST") != 0 )
1538  FeIIPunPop(punch.ipPnunit[ipPun],
1539  true,
1540  (long int)punch.punarg[ipPun][0],
1541  (long int)punch.punarg[ipPun][1],
1542  lgFlag);
1543  }
1544 
1545  /* punch spectra in fits format */
1546  else if( strcmp(punch.chPunch[ipPun],"FITS") == 0 )
1547  {
1548  if( strcmp(chTime,"LAST") == 0 )
1549  {
1551  }
1552  }
1553  /* punch gammas (but without element) */
1554  else if( strcmp(punch.chPunch[ipPun],"GAMt") == 0 )
1555  {
1556  if( strcmp(chTime,"LAST") != 0 )
1557  {
1558  long ns;
1559  /* punch photoionization rates, with the PUNCH GAMMAS command */
1560  for( nelem=0; nelem < LIMELM; nelem++ )
1561  {
1562  for( ion=0; ion <= nelem; ion++ )
1563  {
1564  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1565  {
1566  fprintf( punch.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e",
1567  nelem+1, ion+1, ns+1,
1568  ionbal.PhotoRate_Shell[nelem][ion][ns][0],
1569  ionbal.PhotoRate_Shell[nelem][ion][ns][1] ,
1570  ionbal.PhotoRate_Shell[nelem][ion][ns][2] );
1571 
1572  for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ )
1573  {
1574  fprintf( punch.ipPnunit[ipPun], "%5.2f",
1575  t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) );
1576  }
1577  fprintf( punch.ipPnunit[ipPun], "\n" );
1578  }
1579  }
1580  }
1581  }
1582  }
1583 
1584  /* punch gammas element, ion */
1585  else if( strcmp(punch.chPunch[ipPun],"GAMe") == 0 )
1586  {
1587  if( strcmp(chTime,"LAST") != 0 )
1588  {
1589  int ns;
1590  nelem = (long)punch.punarg[ipPun][0];
1591  ion = (long)punch.punarg[ipPun][1];
1592  /* valence shell */
1593  ns = Heavy.nsShells[nelem][ion]-1;
1594  /* show what some of the ionization sources are */
1595  GammaPrt(
1596  opac.ipElement[nelem][ion][ns][0] ,
1597  opac.ipElement[nelem][ion][ns][1] ,
1598  opac.ipElement[nelem][ion][ns][2] ,
1599  punch.ipPnunit[ipPun],
1600  ionbal.PhotoRate_Shell[nelem][ion][ns][0] ,
1601  ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 );
1602  }
1603  }
1604 
1605  else if( strcmp(punch.chPunch[ipPun],"GAUN") == 0 )
1606  {
1607  /* punch gaunt factors */
1608  if( strcmp(chTime,"LAST") != 0 )
1609  PunchGaunts(punch.ipPnunit[ipPun]);
1610  }
1611 
1612  /* punch parameters of grid calculation */
1613  else if( strcmp(punch.chPunch[ipPun],"GRID") == 0 )
1614  {
1615  /* one time punch of all grid parameters after grid is complete */
1616  if( (strcmp(chTime,"LAST") == 0 ) && grid.lgGridDone )
1617  {
1618  /* start of line gives abort and warning summary */
1619  fprintf(punch.ipPnunit[ipPun], "#Abort?\tWarnings?" );
1620  /* print start of each variable command line */
1621  for( i=0; i< grid.nintparm; i++ )
1622  {
1623  char chStr[10];
1624  strncpy( chStr , optimize.chVarFmt[i] , 9 );
1625  /* make sure this small bit of string is terminated */
1626  chStr[9] = '\0';
1627  fprintf(punch.ipPnunit[ipPun], "\t%s", chStr );
1628  }
1629  fprintf(punch.ipPnunit[ipPun], "\n" );
1630  for( i=0; i<grid.totNumModels; i++ )
1631  {
1632  /* abort / warning summary for this sim */
1633  fprintf(punch.ipPnunit[ipPun], "%c\t%c",
1634  TorF(grid.lgAbort[i]) ,
1635  TorF(grid.lgWarn[i]) );
1636  /* the grid parameters */
1637  for( j=0; j< grid.nintparm; j++ )
1638  {
1639  fprintf(punch.ipPnunit[ipPun], "\t%f", grid.interpParameters[i][j] );
1640  }
1641  fprintf(punch.ipPnunit[ipPun], "\n" );
1642  }
1643  }
1644  }
1645 
1646  else if( strcmp(punch.chPunch[ipPun],"HISp") == 0 )
1647  {
1648  /* punch pressure history of current zone */
1649  if( strcmp(chTime,"LAST") != 0 )
1650  {
1651  /* note if pressure convergence failure occurred in history that follows */
1652  if( !conv.lgConvPres )
1653  {
1654  fprintf( punch.ipPnunit[ipPun],
1655  "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n",
1656  iteration , nzone );
1657  }
1658  /* note if temperature convergence failure occurred in history that follows */
1659  if( !conv.lgConvTemp )
1660  {
1661  fprintf( punch.ipPnunit[ipPun],
1662  "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n",
1663  iteration , nzone );
1664  }
1665  for(i=0; i<conv.hist_pres_npres; ++i )
1666  {
1667  /* punch history of density - pressure, with correct pressure */
1668  fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
1669  iteration,
1670  nzone,
1673  conv.hist_pres_correct[i]);
1674  }
1675  }
1676  }
1677 
1678  else if( strcmp(punch.chPunch[ipPun],"HISt") == 0 )
1679  {
1680  /* punch temperature history of current zone */
1681  if( strcmp(chTime,"LAST") != 0 )
1682  {
1683  /* note if pressure convergence failure occurred in history that follows */
1684  if( !conv.lgConvPres )
1685  {
1686  fprintf( punch.ipPnunit[ipPun],
1687  "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n",
1688  iteration , nzone );
1689  }
1690  /* note if temperature convergence failure occurred in history that follows */
1691  if( !conv.lgConvTemp )
1692  {
1693  fprintf( punch.ipPnunit[ipPun],
1694  "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n",
1695  iteration , nzone );
1696  }
1697  for(i=0; i<conv.hist_temp_ntemp; ++i )
1698  {
1699  /* punch history of density - pressure, with correct pressure */
1700  fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
1701  iteration,
1702  nzone,
1703  conv.hist_temp_temp[i],
1704  conv.hist_temp_heat[i],
1705  conv.hist_temp_cool[i]);
1706  }
1707  }
1708  }
1709 
1710  else if( strncmp(punch.chPunch[ipPun],"H2",2) == 0 )
1711  {
1712  /* all punch info on large H2 molecule include H2 PDR pdr */
1713  H2_PunchDo( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] , chTime, ipPun );
1714  }
1715 
1716  else if( strcmp(punch.chPunch[ipPun],"HEAT") == 0 )
1717  {
1718  /* punch heating, routine in file of same name */
1719  if( strcmp(chTime,"LAST") != 0 )
1720  HeatPunch(punch.ipPnunit[ipPun]);
1721  }
1722 
1723  else if( strncmp(punch.chPunch[ipPun],"HE",2) == 0 )
1724  {
1725  /* various punch helium commands */
1726  /* punch helium line wavelengths */
1727  if( strcmp(punch.chPunch[ipPun] , "HELW") == 0 )
1728  {
1729  if( strcmp(chTime,"LAST") == 0 )
1730  {
1731  /* punch helium & he-like wavelengths, first header */
1732  fprintf( punch.ipPnunit[ipPun],
1733  "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S"
1734  "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" );
1735  fprintf( punch.ipPnunit[ipPun], "\n" );
1736  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
1737  {
1738  /* print element name, nuclear charge */
1739  fprintf( punch.ipPnunit[ipPun], "%li\t%s",
1740  nelem+1 , elementnames.chElementSym[nelem] );
1741  /*prt_wl print floating wavelength in Angstroms, in output format */
1742  fprintf( punch.ipPnunit[ipPun], "\t" );
1743  prt_wl( punch.ipPnunit[ipPun] ,
1744  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].WLAng );
1745  fprintf( punch.ipPnunit[ipPun], "\t" );
1746  prt_wl( punch.ipPnunit[ipPun] ,
1747  Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe1s1S].WLAng );
1748  fprintf( punch.ipPnunit[ipPun], "\t" );
1749  prt_wl( punch.ipPnunit[ipPun] ,
1750  Transitions[ipHE_LIKE][nelem][ipHe2p3P2][ipHe1s1S].WLAng );
1751  fprintf( punch.ipPnunit[ipPun], "\t" );
1752  prt_wl( punch.ipPnunit[ipPun] ,
1753  Transitions[ipHE_LIKE][nelem][ipHe2s3S][ipHe1s1S].WLAng );
1754  fprintf( punch.ipPnunit[ipPun], "\t" );
1755  prt_wl( punch.ipPnunit[ipPun] ,
1756  Transitions[ipHE_LIKE][nelem][ipHe2p3P2][ipHe2s3S].WLAng );
1757  fprintf( punch.ipPnunit[ipPun], "\t" );
1758  prt_wl( punch.ipPnunit[ipPun] ,
1759  Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe2s3S].WLAng );
1760  fprintf( punch.ipPnunit[ipPun], "\t" );
1761  prt_wl( punch.ipPnunit[ipPun] ,
1762  Transitions[ipHE_LIKE][nelem][ipHe2p3P0][ipHe2s3S].WLAng );
1763  fprintf( punch.ipPnunit[ipPun], "\n");
1764  }
1765  }
1766  }
1767  else
1768  TotalInsanity();
1769  }
1770 
1771  /* punch hummer, results needed for Lya transport, to feed into David's routine */
1772  else if( strcmp(punch.chPunch[ipPun],"HUMM") == 0 )
1773  {
1774  eps = Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul/
1775  (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL *dense.eden);
1776  fprintf( punch.ipPnunit[ipPun],
1777  " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1779  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn,
1781  StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*dense.xIonDense[ipHYDROGEN][1],
1782  phycon.te, Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->damp, eps );
1783  }
1784 
1785  else if( strncmp( punch.chPunch[ipPun] , "HYD", 3 ) == 0 )
1786  {
1787  /* various punch hydrogen commands */
1788  if( strcmp(punch.chPunch[ipPun],"HYDc") == 0 )
1789  {
1790  if( strcmp(chTime,"LAST") != 0 )
1791  {
1792  /* punch hydrogen physical conditions */
1793  fprintf( punch.ipPnunit[ipPun],
1794  " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1802  }
1803  }
1804 
1805  else if( strcmp(punch.chPunch[ipPun],"HYDi") == 0 )
1806  {
1807  if( strcmp(chTime,"LAST") != 0 )
1808  {
1809  /* punch hydrogen ionization
1810  * this will be total decays to ground */
1811  RateInter = 0.;
1814  fout = StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop;
1815  /* 06 aug 28, from numLevels_max to _local. */
1816  for( ion=ipH2s; ion < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ion++ )
1817  {
1818  /* this is total decays to ground */
1819  RateInter +=
1820  Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Aul*
1821  (Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pesc +
1822  Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pelec_esc +
1823  Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pdest);
1824  /* total photo from all levels */
1825  fref += iso.gamnc[ipH_LIKE][ipHYDROGEN][ion]*StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop;
1826  /* total col ion from all levels */
1827  stage += iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ion]*dense.eden*
1828  StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop;
1829  fout += StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop;
1830  }
1831  fprintf( punch.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e",
1832  nzone,
1833  /* photo and collision ion rates have units s-1 */
1837 
1838  fprintf( punch.ipPnunit[ipPun], "\t%.2e",
1840 
1841  fprintf( punch.ipPnunit[ipPun],
1842  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1843  dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0],
1844  iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/(ionbal.RateRecomTot[ipHYDROGEN][0]),
1845  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][1][ipRecEsc],
1846  RateInter,
1847  fref/MAX2(1e-37,fout),
1848  stage/MAX2(1e-37,fout),
1849  /* simple H+ */
1850  safe_div( iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][0], dense.eden*dense.xIonDense[ipHYDROGEN][1] ),
1851  secondaries.csupra[ipHYDROGEN][0]);
1852 
1853  GammaPrt(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0],rfield.nflux,iso.ipOpac[ipH_LIKE][ipHYDROGEN][ipH1s],
1854  punch.ipPnunit[ipPun],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*
1855  0.05);
1856  }
1857  }
1858 
1859  else if( strcmp(punch.chPunch[ipPun],"HYDp") == 0 )
1860  {
1861  if( strcmp(chTime,"LAST") != 0 )
1862  {
1863  /* punch hydrogen populations
1864  * first give total atom and ion density [cm-3]*/
1865  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e",
1867  dense.xIonDense[ipHYDROGEN][0],
1868  dense.xIonDense[ipHYDROGEN][1] );
1869 
1870  /* next give state-specific densities [cm-3] */
1871  for( j=ipH1s; j < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]-1; j++ )
1872  {
1873  fprintf( punch.ipPnunit[ipPun], "\t%.2e",
1875  }
1876  fprintf( punch.ipPnunit[ipPun], "\n" );
1877  }
1878  }
1879 
1880  else if( strcmp(punch.chPunch[ipPun],"HYDl") == 0 )
1881  {
1882  if( strcmp(chTime,"LAST") == 0 )
1883  {
1884  /* punch hydrogen line
1885  * gives intensities and optical depths */
1888  for( ipHi=4; ipHi<iso.numLevels_local[ipH_LIKE][ipHYDROGEN]-1; ++ipHi )
1889  {
1890  for( ipLo=ipHi-5; ipLo<ipHi; ++ipLo )
1891  {
1892  if( ipLo<0 )
1893  continue;
1894  fprintf(punch.ipPnunit[ipPun], "%li\t%li\t%.4e\t%.2e\n",
1895  ipHi, ipLo,
1896  /* wl in cm */
1897  1./Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].EnergyWN,
1898  Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Emis->TauIn );
1899  }
1900  }
1901  }
1902  }
1903 
1904  /* punch hydrogen Lya - some details about Lya */
1905  else if( strcmp(punch.chPunch[ipPun],"HYDL") == 0 )
1906  {
1907  if( strcmp(chTime,"LAST") != 0 )
1908  {
1909  /* the population ratio for Lya */
1910  double popul = StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop/SDIV(StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop);
1911  /* the excitation temperature of Lya */
1913  fprintf( punch.ipPnunit[ipPun],
1914  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
1916  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn,
1917  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauTot,
1918  popul,
1919  texc,
1920  phycon.te,
1921  texc/phycon.te ,
1922  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pesc,
1923  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest,
1924  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->pump,
1925  opac.opacity_abs[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1],
1926  opac.albedo[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] );
1927  }
1928  }
1929 
1930  else if( strcmp(punch.chPunch[ipPun],"HYDr") == 0 )
1931  {
1932  /* punch hydrogen recc - recombination cooling for AGN3 */
1933  TempChange(2500.f, false);
1934  while( phycon.te <= 20000. )
1935  {
1936  double r1;
1937  double ThinCoolingCaseB;
1938 
1939  r1 = HydroRecCool(1,0);
1940  ThinCoolingCaseB = pow(10.,((-25.859117 +
1941  0.16229407*phycon.telogn[0] +
1942  0.34912863*phycon.telogn[1] -
1943  0.10615964*phycon.telogn[2])/(1. +
1944  0.050866793*phycon.telogn[0] -
1945  0.014118924*phycon.telogn[1] +
1946  0.0044980897*phycon.telogn[2] +
1947  6.0969594e-5*phycon.telogn[3])))/phycon.te;
1948 
1949  fprintf( punch.ipPnunit[ipPun], " %10.2e\t",
1950  phycon.te);
1951  fprintf( punch.ipPnunit[ipPun], " %10.2e\t",
1952  (r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) );
1953 
1954  fprintf( punch.ipPnunit[ipPun], " %10.2e\t",
1955  r1/(BOLTZMANN*phycon.te));
1956 
1957  fprintf( punch.ipPnunit[ipPun], " %10.2e\n",
1958  ThinCoolingCaseB/(BOLTZMANN*phycon.te));
1959 
1960  TempChange(phycon.te *2.f , false);
1961  }
1962  /* must exit since we have disturbed the solution */
1963  fprintf(ioQQQ , "punch agn now exits since solution is disturbed.\n");
1964  cdEXIT( EXIT_SUCCESS );
1965  }
1966  else
1967  TotalInsanity();
1968  }
1969 
1970  else if( strcmp(punch.chPunch[ipPun],"IONI") == 0 )
1971  {
1972  if( strcmp(chTime,"LAST") == 0 )
1973  {
1974  /* punch mean ionization distribution */
1975  PrtMeanIon( 'i', false , punch.ipPnunit[ipPun] );
1976  }
1977  }
1978 
1979  /* punch ionization rates */
1980  else if( strcmp(punch.chPunch[ipPun],"IONR") == 0 )
1981  {
1982  if( strcmp(chTime,"LAST") != 0 )
1983  {
1984  /* this is element number */
1985  nelem = (long)punch.punarg[ipPun][0];
1986  fprintf( punch.ipPnunit[ipPun],
1987  "%.5e\t%.4e\t%.4e",
1989  dense.eden ,
1990  dynamics.Rate);
1991  /* >>chng 04 oct 15, from nelem+2 to nelem+1 - array over read -
1992  * caught by PnH */
1993  for( ion=0; ion<nelem+1; ++ion )
1994  {
1995  fprintf( punch.ipPnunit[ipPun],
1996  "\t%.4e\t%.4e\t%.4e\t%.4e",
1997  dense.xIonDense[nelem][ion] ,
1998  ionbal.RateIonizTot[nelem][ion] ,
1999  ionbal.RateRecomTot[nelem][ion] ,
2000  dynamics.Source[nelem][ion] );
2001  }
2002  fprintf( punch.ipPnunit[ipPun], "\n");
2003  }
2004  }
2005 
2006  else if( strcmp(punch.chPunch[ipPun]," IP ") == 0 )
2007  {
2008  if( strcmp(chTime,"LAST") == 0 )
2009  {
2010  /* punch valence shell ip's */
2011  for( nelem=0; nelem < LIMELM; nelem++ )
2012  {
2013  int ion_big;
2014  double energy;
2015 
2016  /* this is the largest number of ion stages per line */
2017  const int NELEM_LINE = 10;
2018  /* this loop in case all ions do not fit across page */
2019  for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE )
2020  {
2021  int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem);
2022 
2023  /* new line then element name */
2024  fprintf( punch.ipPnunit[ipPun],
2025  "\n%2.2s", elementnames.chElementSym[nelem]);
2026 
2027  /* print ion stages across line */
2028  for( ion=ion_big; ion <= ion_limit; ++ion )
2029  {
2030  fprintf( punch.ipPnunit[ipPun], "\t%4ld", ion+1 );
2031  }
2032  fprintf( punch.ipPnunit[ipPun], "\n" );
2033 
2034  /* this loop is over all shells */
2035  ASSERT( ion_limit < LIMELM );
2036  /* upper limit is number of shells in atom */
2037  for( ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ )
2038  {
2039 
2040  /* print shell label */
2041  fprintf( punch.ipPnunit[ipPun], "%2.2s", Heavy.chShell[ips]);
2042 
2043  /* loop over possible ions */
2044  for( ion=ion_big; ion<=ion_limit; ++ion )
2045  {
2046 
2047  /* does this subshell exist for this ion? break if it does not*/
2048  /*if( Heavy.nsShells[nelem][ion]<Heavy.nsShells[nelem][0] )*/
2049  if( ips >= Heavy.nsShells[nelem][ion] )
2050  break;
2051 
2052  /* array elements are shell, numb of electrons, element, 0 */
2053  energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0);
2054 
2055  /* now print threshold with correct format */
2056  if( energy < 10. )
2057  {
2058  fprintf( punch.ipPnunit[ipPun], "\t%6.3f", energy );
2059  }
2060  else if( energy < 100. )
2061  {
2062  fprintf( punch.ipPnunit[ipPun], "\t%6.2f", energy );
2063  }
2064  else if( energy < 1000. )
2065  {
2066  fprintf( punch.ipPnunit[ipPun], "\t%6.1f", energy );
2067  }
2068  else
2069  {
2070  fprintf( punch.ipPnunit[ipPun], "\t%6ld", (long)(energy) );
2071  }
2072  }
2073 
2074  /* put cs at end of long line */
2075  fprintf( punch.ipPnunit[ipPun], "\n" );
2076  }
2077  }
2078  }
2079  }
2080  }
2081 
2082  else if( strcmp(punch.chPunch[ipPun],"LINC") == 0 )
2083  {
2084  /* punch line cumulative */
2085  /* lgOK not used, placdholder */
2086  if( strcmp(chTime,"LAST") != 0 )
2087  {
2088  /* not used for anything here, but should not pass baloney*/
2089  lgOK = true;
2090  punch_line(punch.ipPnunit[ipPun],"PUNC",lgOK,chDummy);
2091  }
2092  }
2093 
2094  else if( strcmp(punch.chPunch[ipPun],"LIND") == 0 )
2095  {
2096  /* punch line data, then stop */
2097  PunchLineData(punch.ipPnunit[ipPun]);
2098  }
2099 
2100  else if( strcmp(punch.chPunch[ipPun],"LINL") == 0 )
2101  {
2102  /* punch line labels */
2103  bool lgPrintAll=false;
2104  /* LONG keyword on punch line labels command sets this to 1 */
2105  if( punch.punarg[ipPun][0]>0. )
2106  lgPrintAll = true;
2107  prt_LineLabels(punch.ipPnunit[ipPun] , lgPrintAll );
2108  }
2109 
2110  else if( strcmp(punch.chPunch[ipPun],"LINO") == 0 )
2111  {
2112  if( strcmp(chTime,"LAST") == 0 )
2113  {
2114  /* punch line optical depths, routine is below, file static */
2115  PunchLineStuff(punch.ipPnunit[ipPun],"optical" , punch.punarg[ipPun][0]);
2116  }
2117  }
2118 
2119  else if( strcmp(punch.chPunch[ipPun],"LINP") == 0 )
2120  {
2121  if( strcmp(chTime,"LAST") != 0 )
2122  {
2123  static bool lgFirst=true;
2124  /* punch line populations, need to do this twice if very first
2125  * call since first call to PunchLineStuff generates atomic parameters
2126  * rather than level pops, routine is below, file static */
2127  PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[ipPun][0]);
2128  if( lgFirst )
2129  {
2130  lgFirst = false;
2131  PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[ipPun][0]);
2132  }
2133  }
2134  }
2135 
2136  else if( strcmp(punch.chPunch[ipPun],"LINS") == 0 )
2137  {
2138  /* punch line emissivity */
2139  if( strcmp(chTime,"LAST") != 0 )
2140  {
2141  /* not actually used, but should not pass baloney*/
2142  lgOK = true;
2143  punch_line(punch.ipPnunit[ipPun],"PUNS",lgOK,chDummy);
2144  }
2145  }
2146 
2147  else if( strcmp(punch.chPunch[ipPun],"LINR") == 0 )
2148  {
2149  /* punch line RT */
2150  if( strcmp(chTime,"LAST") != 0 )
2151  {
2152  Punch_Line_RT( punch.ipPnunit[ipPun] , "PUNC" );
2153  }
2154  }
2155 
2156  else if( strcmp(punch.chPunch[ipPun],"LINA") == 0 )
2157  {
2158  /* punch line array */
2159  if( strcmp(chTime,"LAST") == 0 )
2160  {
2161  /* punch out all lines with energies */
2162  for( j=0; j < LineSave.nsum; j++ )
2163  {
2164  if( LineSv[j].wavelength > 0. &&
2165  LineSv[j].sumlin[LineSave.lgLineEmergent] > 0. )
2166  {
2167  /* line energy, in units set with units option */
2168  fprintf( punch.ipPnunit[ipPun], "%12.5e",
2170  /* line label */
2171  fprintf( punch.ipPnunit[ipPun], "\t%4.4s ",
2172  LineSv[j].chALab );
2173  /* wavelength */
2174  prt_wl( punch.ipPnunit[ipPun], LineSv[j].wavelength );
2175  /* intrinsic intensity */
2176  fprintf( punch.ipPnunit[ipPun], "\t%8.3f",
2177  log10(SDIV(LineSv[j].sumlin[0]) ) + radius.Conv2PrtInten );
2178  /* emergent line intensity, r recombination */
2179  fprintf( punch.ipPnunit[ipPun], "\t%8.3f",
2180  log10(SDIV(LineSv[j].sumlin[1]) ) + radius.Conv2PrtInten );
2181  /* type of line, i for info, etc */
2182  fprintf( punch.ipPnunit[ipPun], " \t%c\n",
2183  LineSv[j].chSumTyp);
2184  }
2185  }
2186  }
2187  }
2188 
2189  else if( strcmp(punch.chPunch[ipPun],"LINI") == 0 )
2190  {
2191  if( strcmp(chTime,"LAST") == 0 &&
2193  {
2194  /* this is the last zone
2195  * punch line intensities - but do not do last zone twice */
2197  }
2198  else if( strcmp(chTime,"LAST") != 0 )
2199  {
2200  /* following so we only punch first zone if LinEvery reset */
2201  if( (punch.lgLinEvery && nzone == 1) ||
2203  nzone )
2204  {
2205  /* this is middle of calculation
2206  * punch line intensities */
2208  }
2209  }
2210  }
2211 
2212  else if( strcmp( punch.chPunch[ipPun],"LEIL") == 0)
2213  {
2214  /* some line intensities for the Leiden PDR,
2215  * but only do this when calculation is complete */
2216  if( strcmp(chTime,"LAST") == 0 )
2217  {
2218  double absval , rel;
2219  long int n;
2220  /* the lines we will find,
2221  * for a sample list of PDR lines look at LineList_PDR_H2.dat
2222  * in the cloudy data dir */
2223  /* the number of H2 lines */
2224  const int NLINE_H2 = 31;
2225  /* the number of lines which are not H2 */
2226  const int NLINE_NOTH_H2 = 5;
2227  /* the labels and wavelengths for the lines that are not H2 */
2228  char chLabel[NLINE_NOTH_H2][5]=
2229  {"C 2", "O 1", "O 1","C 1", "C 1" };
2230  double Wl[NLINE_NOTH_H2]=
2231  {157.6 , 63.17 , 145.5 ,609.2 , 369.7 , };
2232  /* these are wavelengths in microns, conv to Angstroms before call */
2233  /* >>chng 05 sep 06, many of following wavelengths updated to agree
2234  * with output - apparently not updated when energies changed */
2235  double Wl_H2[NLINE_H2]=
2236  {2.121 ,
2237  28.21 , 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 ,
2238  4.408 , 4.180 , 3.996 , 3.845 , 3.724 , 3.626 , 3.547 , 3.485 , 3.437 , 3.403 ,
2239  3.381 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603};
2240  /* print a header for the lines */
2241  for( n=0; n<NLINE_NOTH_H2; ++n )
2242  {
2243  fprintf(punch.ipPnunit[ipPun], "%s\t%.2f",chLabel[n] , Wl[n]);
2244  /* get the line, non positive return says didn't find it */
2245  /* arguments are 4-char label, wavelength, return log total intensity, linear rel inten */
2246  if( cdLine( chLabel[n] , (realnum)(Wl[n]*1e4) , &absval , &rel ) <= 0 )
2247  {
2248  fprintf(punch.ipPnunit[ipPun], " did not find\n");
2249  }
2250  else
2251  {
2252  fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
2253  }
2254  }
2255  fprintf(punch.ipPnunit[ipPun], "\n\n\n");
2256 
2257  /* only print the H2 lines if the big molecule is turned on */
2258  if( h2.lgH2ON )
2259  {
2260  fprintf(punch.ipPnunit[ipPun],
2261  "Here are some of the H2 Intensities, The first one is the\n"
2262  "1-0 S(0) line and the following ones are the 0-0 S(X)\n"
2263  "lines where X goes from 0 to 29\n\n");
2264  for( n=0; n<NLINE_H2; ++n )
2265  {
2266  fprintf(punch.ipPnunit[ipPun], "%s\t%.3f","H2 " , Wl_H2[n]);
2267  /* get the line, non positive return says didn't find it */
2268  if( cdLine( "H2 " , (realnum)(Wl_H2[n]*1e4) , &absval , &rel ) <= 0 )
2269  {
2270  fprintf(punch.ipPnunit[ipPun], " did not find\n");
2271  }
2272  else
2273  {
2274  fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
2275  }
2276  }
2277  }
2278  }
2279  }
2280 
2281  else if( strcmp( punch.chPunch[ipPun],"LEIS") == 0)
2282  {
2283  if( strcmp(chTime,"LAST") != 0 )
2284  {
2285  /* get some column densities we shall need */
2286  double col_ci , col_oi , col_cii, col_heii;
2287  if( cdColm("carb" , 1 , &col_ci ) )
2288  TotalInsanity();
2289  if( cdColm("carb" , 2 , &col_cii ) )
2290  TotalInsanity();
2291  if( cdColm("oxyg" , 1 , &col_oi ) )
2292  TotalInsanity();
2293  if( cdColm("heli" , 2 , &col_heii ) )
2294  TotalInsanity();
2295  /* punch Leiden structure - some numbers for the Leiden PDR model comparisons */
2296  fprintf( punch.ipPnunit[ipPun],
2297  "%.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"
2298  "%.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"
2299  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
2300  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
2301  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2302  /* depth of this point */
2304  /* A_V for an extended source */
2305  0.00,
2306  /* A_V for a point source */
2308  /* temperature */
2309  phycon.te ,
2311  hmi.H2_total,
2312  dense.xIonDense[ipCARBON][0],
2313  dense.xIonDense[ipCARBON][1],
2314  dense.xIonDense[ipOXYGEN][0],
2315  findspecies("CO")->hevmol,
2316  findspecies("O2")->hevmol,
2317  findspecies("CH")->hevmol,
2318  findspecies("OH")->hevmol,
2319  dense.eden,
2320  dense.xIonDense[ipHELIUM][1],
2322  hmi.Hmolec[ipMH3p],
2325  col_ci,
2326  col_cii,
2327  col_oi,
2328  findspecies("CO")->hevcol,
2329  findspecies("O2")->hevcol,
2330  findspecies("CH")->hevcol,
2331  findspecies("OH")->hevcol,
2333  col_heii,
2339  /* CO and C dissociation rate */
2340  CO_findrk("PHOTON,CO=>C,O"),
2341  /* total CI ionization rate */
2342  ionbal.PhotoRate_Shell[ipCARBON][0][2][0],
2343  /* total heating, erg cm-3 s-1 */
2344  thermal.htot,
2345  /* total cooling, erg cm-3 s-1 */
2346  thermal.ctot,
2347  /* GrnP grain photo heating */
2348  thermal.heating[0][13],
2349  /* grain collisional cooling */
2350  MAX2(0.,gv.GasCoolColl),
2351  /* grain collisional heating */
2352  -1.*MIN2(0.,gv.GasCoolColl),
2353  /* COds - CO dissociation heating */
2354  thermal.heating[0][9],
2355  /* H2dH-Heating due to H2 dissociation */
2357  /* H2vH-Heating due to collisions with H2 */
2359  /* ChaT - charge transfer heating */
2360  thermal.heating[0][24] ,
2361  /* cosmic ray heating */
2362  thermal.heating[1][6] ,
2363  /* heating due to atoms of various heavy elements */
2367  thermal.heating[ipIRON][0],
2368  thermal.heating[ipSODIUM][0],
2370  thermal.heating[ipCARBON][0],
2375  TauLines[ipT146].Coll.cool );
2376  }
2377  }
2378 
2379  else if( strcmp( punch.chPunch[ipPun],"LLST") == 0)
2380  {
2381  /* punch linelist command - do on last iteration */
2382  if( strcmp(chTime,"LAST") == 0 )
2383  {
2384  fprintf( punch.ipPnunit[ipPun],
2385  "%li" , iteration );
2386 
2387  /* -1 is flag saying that this punch command was not set */
2388  if( punch.nLineList[ipPun] < 0 )
2389  TotalInsanity();
2390 
2391  /* loop over all lines in the file we read */
2392  for( j=0; j<punch.nLineList[ipPun]; ++j )
2393  {
2394  double relative , absolute, PrtQuantity;
2395  if( (cdLine( punch.chLineListLabel[ipPun][j] ,
2396  punch.wlLineList[ipPun][j] ,
2397  &relative , &absolute ))<=0 )
2398  {
2399  if( !h2.lgH2ON && strncmp( punch.chLineListLabel[ipPun][j] , "H2 " , 4 )==0 )
2400  {
2401  static bool lgMustPrintFirstTime = true;
2402  if( lgMustPrintFirstTime )
2403  {
2404  /* it's an H2 line and H2 is not being done - ignore it */
2405  fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
2406  "included, so I will ignore it. Log intensity set to -30.\n" );
2407  fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
2408  lgMustPrintFirstTime = false;
2409  }
2410  relative = -30.f;
2411  absolute = -30.f;
2412  }
2413  else if( lgAbort )
2414  {
2415  /* we are in abort mode */
2416  relative = -30.f;
2417  absolute = -30.f;
2418  }
2419  else
2420  {
2421  fprintf(ioQQQ,"DISASTER - did not find a line in the Line List table\n");
2422  cdEXIT(EXIT_FAILURE);
2423  }
2424  }
2425 
2426  /* options to do either relative or absolute intensity
2427  * default is relative, is absolute keyword on line then
2428  * punarg set to 1 */
2429  /* straight line intensities */
2430  if( punch.punarg[ipPun][0] > 0 )
2431  PrtQuantity = absolute;
2432  else
2433  PrtQuantity = relative;
2434  /* if taking ratio print every other line as ratio
2435  * with previous line */
2436  if( punch.lgLineListRatio[ipPun] )
2437  {
2438  /* do line pair ratios */
2439  static double SaveQuantity = 0;
2440  if( is_odd(j) )
2441  fprintf( punch.ipPnunit[ipPun], "\t%.3e" ,
2442  SaveQuantity / SDIV( PrtQuantity ) );
2443  else
2444  SaveQuantity = PrtQuantity;
2445  }
2446  else
2447  {
2448  fprintf( punch.ipPnunit[ipPun], "\t%.3e" , PrtQuantity );
2449  }
2450  }
2451  fprintf( punch.ipPnunit[ipPun], "\n" );
2452  }
2453  }
2454 
2455  else if( strcmp( punch.chPunch[ipPun],"RCO ") == 0)
2456  {
2457  /* punch chemistry CO rates command */
2458  if( strcmp(chTime,"LAST") != 0 )
2459  {
2460  CO_punch_mol(punch.ipPnunit[ipPun],"CO",NULL,radius.depth_mid_zone);
2461  }
2462  }
2463 
2464  /*>>chng 06 jul 13, Nick Abel add OH rates */
2465  else if( strcmp( punch.chPunch[ipPun],"ROH ") == 0)
2466  {
2467  /* punch chemistry OH rates command */
2468  if( strcmp(chTime,"LAST") != 0 )
2469  {
2470  CO_punch_mol(punch.ipPnunit[ipPun],"OH",NULL,radius.depth_mid_zone);
2471  }
2472  }
2473  else if( strcmp(punch.chPunch[ipPun],"MAP ") == 0 )
2474  {
2475  /* do the map now if we are at the zone, or if this
2476  * is the LAST call to this routine and map not done yet */
2477  if( !hcmap.lgMapDone &&
2478  (nzone == hcmap.MapZone || strcmp(chTime,"LAST") == 0 ) )
2479  {
2480  lgTlkSav = called.lgTalk;
2481  called.lgTalk = true;
2482  hcmap.lgMapBeingDone = true;
2483  map_do(punch.ipPnunit[ipPun] , " map");
2484  called.lgTalk = lgTlkSav;
2485  }
2486  }
2487 
2488  else if( strcmp(punch.chPunch[ipPun],"MOLE") == 0 )
2489  {
2490  if( strcmp(chTime,"LAST") != 0 )
2491  {
2492  static bool lgMustPrintHeader = true;
2493  if( lgMustPrintHeader )
2494  {
2495  lgMustPrintHeader = false;
2496  fprintf( punch.ipPnunit[ipPun],
2497  "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate");
2498 
2499  for(i=0; i<N_H_MOLEC; ++i )
2500  {
2501  fprintf( punch.ipPnunit[ipPun], "\t%s", hmi.chLab[i] );
2502  }
2503  for(i=0; i<mole.num_comole_calc; ++i )
2504  {
2505  fprintf( punch.ipPnunit[ipPun], "\t%s",COmole[i]->label );
2506  }
2507  fprintf( punch.ipPnunit[ipPun], "\n");
2508  }
2509  /* molecules, especially for PDR, first give radius */
2510  fprintf( punch.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone );
2511 
2512  /* visual extinction of point source (star)*/
2513  fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
2514 
2515  /* visual extinction of an extended source (like a PDR)*/
2516  fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
2517 
2518  /* carbon monoxide photodissociation rate */
2519  fprintf( punch.ipPnunit[ipPun], "%.5e\t" , CO_findrk("PHOTON,CO=>C,O") );
2520 
2521  /* carbon recombination rate */
2522  fprintf( punch.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] );
2523 
2524  /* now do all the H molecules */
2525  for(j=0; j<N_H_MOLEC; ++j )
2526  {
2527  fprintf(punch.ipPnunit[ipPun],"\t%.2e",hmi.Hmolec[j] );
2528  }
2529 
2530  /* now all the C molecules */
2531  for(j=0; j<mole.num_comole_calc; ++j )
2532  {
2533  fprintf(punch.ipPnunit[ipPun],"\t%.2e",COmole[j]->hevmol );
2534  }
2535  fprintf(punch.ipPnunit[ipPun],"\n");
2536  }
2537  }
2538 
2539  else if( strcmp(punch.chPunch[ipPun],"OPAC") == 0 )
2540  {
2541  /* punch opacity - routine will parse which type of opacity punch to do */
2542  if( strcmp(chTime,"LAST") == 0 )
2543  punch_opacity(punch.ipPnunit[ipPun],ipPun);
2544  }
2545 
2546  /* punch coarse optical depths command */
2547  else if( strcmp(punch.chPunch[ipPun],"OPTc") == 0 )
2548  {
2549  if( strcmp(chTime,"LAST") == 0 )
2550  {
2551  for( j=0; j < rfield.nflux; j++ )
2552  {
2553  fprintf( punch.ipPnunit[ipPun],
2554  "%12.4e\t%.3e\t%12.4e\t%.3e\n",
2555  AnuUnit(rfield.AnuOrg[j]),
2557  opac.TauAbsFace[j],
2558  opac.TauScatFace[j] );
2559  }
2560  }
2561  }
2562 
2563  /* punch fine optical depths command */
2564  else if( strcmp(punch.chPunch[ipPun],"OPTf") == 0 )
2565  {
2566  if( strcmp(chTime,"LAST") == 0 )
2567  {
2568  long nu_hi , nskip;
2569  if( punch.punarg[ipPun][0] > 0. )
2570  /* possible lower bounds to energy range - will be zero if not set */
2571  j = ipFineCont( punch.punarg[ipPun][0] );
2572  else
2573  j = 0;
2574 
2575  /* upper limit */
2576  if( punch.punarg[ipPun][1]> 0. )
2577  nu_hi = ipFineCont( punch.punarg[ipPun][1]);
2578  else
2579  nu_hi = rfield.nfine;
2580 
2581  /* we will bring nskip cells together into one printed
2582  * number to make output smaller - default is 10 */
2583  nskip = (long)punch.punarg[ipPun][2];
2584  do
2585  {
2586  realnum sum1 = rfield.fine_opt_depth[j];
2587  realnum sum2 = rfield.fine_opac_zone[j];
2588  /* want to report the central wavelength of the cell */
2589  realnum xnu = rfield.fine_anu[j];
2590  for( jj=1; jj<nskip; ++jj )
2591  {
2592  sum1 += rfield.fine_opt_depth[j+jj];
2593  sum2 += rfield.fine_opac_zone[j+jj];
2594  xnu += rfield.fine_anu[j+jj];
2595  }
2596  if( sum2 > 0. )
2597  fprintf( punch.ipPnunit[ipPun],
2598  "%12.6e\t%.3e\t%.3e\n",
2599  AnuUnit(xnu/nskip),
2600  sum1/nskip ,
2601  sum2/nskip);
2602  j += nskip;
2603  }while( j < nu_hi );
2604  }
2605  }
2606 
2607  else if( strcmp(punch.chPunch[ipPun]," OTS") == 0 )
2608  {
2609  ConMax = 0.;
2610  xLinMax = 0.;
2611  opConSum = 0.;
2612  opLinSum = 0.;
2613  ipLinMax = 1;
2614  ipConMax = 1;
2615 
2616  for( j=0; j < rfield.nflux; j++ )
2617  {
2618  opConSum += rfield.otscon[j]*opac.opacity_abs[j];
2619  opLinSum += rfield.otslin[j]*opac.opacity_abs[j];
2620  if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax )
2621  {
2622  xLinMax = rfield.otslin[j]*opac.opacity_abs[j];
2623  ipLinMax = j+1;
2624  }
2625  if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax )
2626  {
2627  ConMax = rfield.otscon[j]*opac.opacity_abs[j];
2628  ipConMax = j+1;
2629  }
2630  }
2631  fprintf( punch.ipPnunit[ipPun],
2632  "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n",
2633  opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1]
2634  , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1]
2635  , rfield.anu[ipConMax-1], ConMax );
2636  }
2637 
2638  else if( strcmp(punch.chPunch[ipPun],"OVER") == 0 )
2639  {
2640  /* punch overview
2641  * this is the floor for the smallest ionization fractions printed */
2642  double toosmall = SMALLFLOAT ,
2643  hold;
2644 
2645  /* overview of model results,
2646  * depth, te, hden, eden, ion fractions H, He, c, O */
2647  if( strcmp(chTime,"LAST") != 0 )
2648  {
2649 
2650  /* print the depth */
2651  fprintf( punch.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone );
2652 
2653  /* temperature, heating */
2654  if(dynamics.Cool > dynamics.Heat)
2655  {
2656  fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f",
2657  log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat)) );
2658  }
2659  else
2660  {
2661  double diff = fabs(thermal.htot-dynamics.Cool);
2662  fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f",
2663  log10(phycon.te), log10(SDIV(diff)) );
2664  }
2665 
2666  /* hydrogen and electron densities */
2667  fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f",
2668  log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) );
2669 
2670  /* molecular fraction of hydrogen */
2671  fprintf( punch.ipPnunit[ipPun], "\t%.4f",
2672  /*log10(MAX2(toosmall,2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN])) );*/
2673  log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) );
2674 
2675  /* ionization fractions of hydrogen */
2676  fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f",
2677  log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])),
2678  log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );
2679 
2680  /* ionization fractions of helium */
2681  for( j=1; j <= 3; j++ )
2682  {
2683  double arg1 = SDIV(dense.gas_phase[ipHELIUM]);
2684  arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 );
2685  fprintf( punch.ipPnunit[ipPun], "\t%.4f",
2686  log10(arg1) );
2687  }
2688 
2689  /* carbon monoxide molecular fraction of CO */
2690  hold = SDIV(dense.gas_phase[ipCARBON]);
2691  hold = findspecies("CO")->hevmol/hold;
2692  hold = MAX2(toosmall, hold );
2693  fprintf( punch.ipPnunit[ipPun], "\t%.4f", log10(hold) );
2694 
2695  /* ionization fractions of carbon */
2696  for( j=1; j <= 4; j++ )
2697  {
2698  hold = SDIV(dense.gas_phase[ipCARBON]);
2699  hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold);
2700  fprintf( punch.ipPnunit[ipPun], "\t%.4f",
2701  log10(hold) );
2702  }
2703 
2704  /* ionization fractions of oxygen */
2705  for( j=1; j <= 6; j++ )
2706  {
2707  hold = SDIV(dense.gas_phase[ipOXYGEN]);
2708  hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold);
2709  fprintf( punch.ipPnunit[ipPun], "\t%.4f",
2710  log10(hold) );
2711  }
2712 
2713  /* visual extinction of point source (star)*/
2714  fprintf( punch.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point);
2715 
2716  /* visual extinction of an extended source (like a PDR)*/
2717  fprintf( punch.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended);
2718  }
2719  }
2720 
2721  else if( strcmp(punch.chPunch[ipPun]," PDR") == 0 )
2722  {
2723  /* this is the punch PDR command */
2724  if( strcmp(chTime,"LAST") != 0 )
2725  {
2726  /* convert optical depth at wavelength of V filter
2727  * into magnitudes of extinction */
2728  /* >>chyng 03 feb 25, report extinction to illuminated face,
2729  * rather than total extinction which included far side when
2730  * sphere was set */
2731  /*av = opac.TauTotalGeo[0][rfield.ipV_filter-1]*1.08574;*/
2732 
2733  fprintf( punch.ipPnunit[ipPun],
2734  "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t",
2736  /* total hydrogen column density, all forms */
2738  phycon.te,
2739  /* fraction of H that is atomic */
2741  /* ratio of n(H2) to total H, == 0.5 when fully molecular */
2744  /* atomic to total carbon */
2746  findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
2747  findspecies("H2O")->hevmol/SDIV(dense.gas_phase[ipOXYGEN]),
2748  /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */
2750 
2751  /* visual extinction due to dust alone, of point source (star)*/
2752  fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
2753 
2754  /* visual extinction due to dust alone, of an extended source (like a PDR)*/
2755  fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
2756 
2757  /* visual extinction (all sources) of a point source (like a PDR)*/
2758  fprintf( punch.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] );
2759  }
2760  }
2761 
2762  else if( strcmp(punch.chPunch[ipPun],"PHYS") == 0 )
2763  {
2764  if( strcmp(chTime,"LAST") != 0 )
2765  {
2766  /* punch physical conditions */
2767  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2770  }
2771  }
2772 
2773  else if( strcmp(punch.chPunch[ipPun],"PRES") == 0 )
2774  {
2775  /* the punch pressure command */
2776  if( strcmp(chTime,"LAST") != 0 )
2777  {
2778  fprintf( punch.ipPnunit[ipPun],
2779  "%.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",
2780  /*A 1 #P depth */
2782  /*B 2 Pcorrect */
2784  /*C 3 Pcurrent */
2786  /*D 4 Pln + pintg
2787  * >>chng 06 apr 19, subtract pinzon the acceleration added in this zone
2788  * since is not total at outer edge of zone, above is at inner edge */
2790  /*E 5 pgas (0) */
2792  /*F 6 Pgas */
2794  /*G 7 Pram */
2796  /*H 8 P rad in lines */
2798  /*I 9 Pinteg subtract continuum rad pres which has already been added on */
2800  /*J 10 V(wind km/s) wind speed in km/s */
2801  wind.windv/1e5,
2802  /*K cad(km/s) sound speed in km/s */
2804  /* the magnetic pressure */
2805  magnetic.pressure ,
2806  /* the local turbulent velocity in km/s */
2807  DoppVel.TurbVel/1e5 ,
2808  /* turbulent pressure */
2810  TorF(conv.lgConvPres) );
2811  }
2812  }
2813 
2814  else if( punch.chPunch[ipPun][0]=='R' )
2815  {
2816  /* work around internal limits to Microsoft vs compiler */
2817  if( strcmp(punch.chPunch[ipPun],"RADI") == 0 )
2818  {
2819  /* punch radius information for all zones */
2820  if( strcmp(chTime,"LAST") != 0 )
2821  {
2822  fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
2824  radius.drad );
2825  }
2826  }
2827 
2828  else if( strcmp(punch.chPunch[ipPun],"RADO") == 0 )
2829  {
2830  /* punch radius information for only the last zone */
2831  if( strcmp(chTime,"LAST") == 0 )
2832  {
2833  fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
2835  radius.drad );
2836  }
2837  }
2838 
2839  else if( strcmp(punch.chPunch[ipPun],"RESU") == 0 )
2840  {
2841  /* punch results of the calculation */
2842  if( strcmp(chTime,"LAST") == 0 )
2843  punResults(punch.ipPnunit[ipPun]);
2844  }
2845  else
2846  {
2847  /* this can't happen */
2848  TotalInsanity();
2849  }
2850  }
2851 
2852  else if( strcmp(punch.chPunch[ipPun],"SECO") == 0 )
2853  {
2854  /* punch secondary ionization */
2855  if( strcmp(chTime,"LAST") != 0 )
2856  fprintf(punch.ipPnunit[ipPun],
2857  "%.5e\t%.3e\t%.3e\t%.3e\n",
2858  radius.depth ,
2860  secondaries.csupra[ipHYDROGEN][0]*2.02,
2861  secondaries.x12tot );
2862  }
2863 
2864  else if( strcmp(punch.chPunch[ipPun],"SOUS") == 0 )
2865  {
2866  /* full spectrum of continuum source function at 1 depth
2867  * command was "punch source spectrum" */
2868  if( strcmp(chTime,"LAST") != 0 )
2869  {
2870  limit = MIN2(rfield.ipMaxBolt,rfield.nflux);
2871  for( j=0; j < limit; j++ )
2872  {
2873  fprintf( punch.ipPnunit[ipPun],
2874  "%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
2875  rfield.anu[j],
2877  opac.opacity_abs[j],
2880  }
2881  }
2882  }
2883 
2884  else if( strcmp(punch.chPunch[ipPun],"SOUD") == 0 )
2885  {
2886  /* parts of continuum source function vs depth
2887  * command was punch source function depth */
2889  fprintf( punch.ipPnunit[ipPun],
2890  "%.4e\t%.4e\t%.4e\t%.4e\n",
2891  opac.TauAbsFace[j-1],
2892  rfield.ConEmitLocal[nzone][j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]),
2895  }
2896 
2897  /* this is punch special option */
2898  else if( strcmp(punch.chPunch[ipPun],"SPEC") == 0 )
2899  {
2900  PunchSpecial(punch.ipPnunit[ipPun],chTime);
2901  }
2902 
2903  else if( strcmp(punch.chPunch[ipPun],"TEGR") == 0 )
2904  {
2905  /* punch history of heating and cooling */
2906  if( strcmp(chTime,"LAST") != 0 )
2907  {
2908  fprintf( punch.ipPnunit[ipPun], " " );
2909  for( j=0; j < NGRID; j++ )
2910  {
2911  fprintf( punch.ipPnunit[ipPun],
2912  "%4ld%11.3e%11.3e%11.3e\n",
2913  thermal.nZonGrid[j],
2914  thermal.TeGrid[j],
2915  thermal.HtGrid[j],
2916  thermal.ClGrid[j] );
2917  }
2918  }
2919  }
2920 
2921  else if( strcmp(punch.chPunch[ipPun],"TEMP") == 0 )
2922  {
2923  static double deriv_old=-1;
2924  double deriv=-1. , deriv_sec;
2925  /* temperature and its derivatives */
2926  fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e",
2928  phycon.te,
2929  thermal.dCooldT );
2930  /* if second zone then have one deriv */
2931  if( nzone >1 )
2932  {
2933  deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
2934  fprintf( punch.ipPnunit[ipPun], "\t%.2e", deriv );
2935  /* if third zone then have second deriv */
2936  if( nzone > 2 )
2937  {
2938  deriv_sec = (deriv-deriv_old)/ radius.drad;
2939  fprintf( punch.ipPnunit[ipPun], "\t%.2e",
2940  deriv_sec );
2941  }
2942  deriv_old = deriv;
2943  }
2944  fprintf( punch.ipPnunit[ipPun], "\n");
2945  }
2946 
2947  /* time dependent model */
2948  else if( strcmp(punch.chPunch[ipPun],"TIMD") == 0 )
2949  {
2950  if( strcmp(chTime,"LAST") == 0 )
2951  DynaPunchTimeDep( punch.ipPnunit[ipPun] , "END" );
2952  }
2953 
2954  /* execution time per zone */
2955  else if( strcmp(punch.chPunch[ipPun],"XTIM") == 0 )
2956  {
2957  static double ElapsedTime , ZoneTime;
2958  if( nzone<=1 )
2959  {
2960  ElapsedTime = cdExecTime();
2961  ZoneTime = 0.;
2962  }
2963  else
2964  {
2965  double t = cdExecTime();
2966  ZoneTime = t - ElapsedTime;
2967  ElapsedTime = t;
2968  }
2969 
2970  /* zone, time for this zone, elapsed time */
2971  fprintf( punch.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n",
2972  nzone, ZoneTime , ElapsedTime );
2973  }
2974 
2975  else if( strcmp(punch.chPunch[ipPun],"TPRE") == 0 )
2976  {
2977  /* temperature and its predictors, turned on with punch tprid */
2978  fprintf( punch.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n",
2980  (phycon.TeProp- phycon.te)/phycon.te );
2981  }
2982 
2983  else if( strcmp(punch.chPunch[ipPun],"WIND") == 0 )
2984  {
2985  /* wind velocity, radiative acceleration, and ratio total
2986  * to electron scattering acceleration */
2987  /* first test only punch last zone */
2988  if( (punch.punarg[ipPun][0] == 0 && strcmp(chTime,"LAST") == 0)
2989  ||
2990  /* this test punch all zones */
2991  (punch.punarg[ipPun][0] == 1 && strcmp(chTime,"LAST") != 0 ) )
2992  {
2993  fprintf( punch.ipPnunit[ipPun],
2994  "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
2997  wind.windv,
2998  wind.AccelTot,
2999  wind.AccelLine,
3000  wind.AccelCont ,
3001  wind.fmul ,
3002  wind.AccelGravity );
3003  }
3004  }
3005 
3006  else if( strcmp(punch.chPunch[ipPun],"XATT") == 0 )
3007  {
3008  /* attenuated incident continuum */
3009  ASSERT( grid.lgOutputTypeOn[2] );
3010 
3011  if( strcmp(chTime,"LAST") == 0 )
3012  {
3013  if( grid.lgGrid )
3014  punchFITSfile( punch.ipPnunit[ipPun], 2 );
3015  else
3016  {
3017  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3018  cdEXIT(EXIT_FAILURE);
3019  }
3020  }
3021  }
3022  else if( strcmp(punch.chPunch[ipPun],"XRFI") == 0 )
3023  {
3024  /* reflected incident continuum */
3025  ASSERT( grid.lgOutputTypeOn[3] );
3026 
3027  if( strcmp(chTime,"LAST") == 0 )
3028  {
3029  if( grid.lgGrid )
3030  punchFITSfile( punch.ipPnunit[ipPun], 3 );
3031  else
3032  {
3033  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3034  cdEXIT(EXIT_FAILURE);
3035  }
3036  }
3037  }
3038  else if( strcmp(punch.chPunch[ipPun],"XINC") == 0 )
3039  {
3040  /* incident continuum */
3041  ASSERT( grid.lgOutputTypeOn[1] );
3042 
3043  if( strcmp(chTime,"LAST") == 0 )
3044  {
3045  if( grid.lgGrid )
3046  punchFITSfile( punch.ipPnunit[ipPun], 1 );
3047  else
3048  {
3049  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3050  cdEXIT(EXIT_FAILURE);
3051  }
3052  }
3053  }
3054  else if( strcmp(punch.chPunch[ipPun],"XDFR") == 0 )
3055  {
3056  /* reflected diffuse continuous emission */
3057  ASSERT( grid.lgOutputTypeOn[5] );
3058 
3059  if( strcmp(chTime,"LAST") == 0 )
3060  {
3061  if( grid.lgGrid )
3062  punchFITSfile( punch.ipPnunit[ipPun], 5 );
3063  else
3064  {
3065  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3066  cdEXIT(EXIT_FAILURE);
3067  }
3068  }
3069  }
3070  else if( strcmp(punch.chPunch[ipPun],"XDFO") == 0 )
3071  {
3072  /* diffuse continuous emission outward */
3073  ASSERT( grid.lgOutputTypeOn[4] );
3074 
3075  if( strcmp(chTime,"LAST") == 0 )
3076  {
3077  if( grid.lgGrid )
3078  punchFITSfile( punch.ipPnunit[ipPun], 4 );
3079  else
3080  {
3081  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3082  cdEXIT(EXIT_FAILURE);
3083  }
3084  }
3085  }
3086  else if( strcmp(punch.chPunch[ipPun],"XLNR") == 0 )
3087  {
3088  /* reflected lines */
3089  ASSERT( grid.lgOutputTypeOn[7] );
3090 
3091  if( strcmp(chTime,"LAST") == 0 )
3092  {
3093  if( grid.lgGrid )
3094  punchFITSfile( punch.ipPnunit[ipPun], 7 );
3095  else
3096  {
3097  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3098  cdEXIT(EXIT_FAILURE);
3099  }
3100  }
3101  }
3102  else if( strcmp(punch.chPunch[ipPun],"XLNO") == 0 )
3103  {
3104  /* outward lines */
3105  ASSERT( grid.lgOutputTypeOn[6] );
3106 
3107  if( strcmp(chTime,"LAST") == 0 )
3108  {
3109  if( grid.lgGrid )
3110  punchFITSfile( punch.ipPnunit[ipPun], 6 );
3111  else
3112  {
3113  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3114  cdEXIT(EXIT_FAILURE);
3115  }
3116  }
3117  }
3118  else if( strcmp(punch.chPunch[ipPun],"XREF") == 0 )
3119  {
3120  /* total reflected, lines and continuum */
3121  ASSERT( grid.lgOutputTypeOn[9] );
3122 
3123  if( strcmp(chTime,"LAST") == 0 )
3124  {
3125  if( grid.lgGrid )
3126  punchFITSfile( punch.ipPnunit[ipPun], 9 );
3127  else
3128  {
3129  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3130  cdEXIT(EXIT_FAILURE);
3131  }
3132  }
3133  }
3134  else if( strcmp(punch.chPunch[ipPun],"XTRN") == 0 )
3135  {
3136  /* total outward, lines and continuum */
3137  ASSERT( grid.lgOutputTypeOn[8] );
3138 
3139  if( strcmp(chTime,"LAST") == 0 )
3140  {
3141  if( grid.lgGrid )
3142  punchFITSfile( punch.ipPnunit[ipPun], 8 );
3143  else
3144  {
3145  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3146  cdEXIT(EXIT_FAILURE);
3147  }
3148  }
3149  }
3150  else if( strcmp(punch.chPunch[ipPun],"XSPM") == 0 )
3151  {
3152  /* exp(-tau) to the illuminated face */
3153  ASSERT( grid.lgOutputTypeOn[10] );
3154 
3155  if( strcmp(chTime,"LAST") == 0 )
3156  {
3157  if( grid.lgGrid )
3158  punchFITSfile( punch.ipPnunit[ipPun], 10 );
3159  else
3160  {
3161  fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
3162  cdEXIT(EXIT_FAILURE);
3163  }
3164  }
3165  }
3166  /* there are a few "punch" commands that are handled elsewhere
3167  * punch dr is an example. These will have lgReadPunch set
3168  * false */
3169  else if( punch.lgRealPunch[ipPun] )
3170  {
3171  /* this is insanity, internal flag set in ParsePunch not seen here */
3172  fprintf( ioQQQ, " PROBLEM DISASTER PunchDo does not recognize flag %4.4s set by ParsePunch. This is impossible.\n",
3173  punch.chPunch[ipPun] );
3174  TotalInsanity();
3175  }
3176 
3177  /* print special hash string to separate out various iterations
3178  * chTime is LAST on last iteration
3179  * punch.lgHashEndIter flag is true by default, set false
3180  * with "no hash" keyword on punch command
3181  * punch.lg_separate_iterations is true by default, set false
3182  * when punch time dependent calcs since do not want special
3183  * character between time steps
3184  * grid.lgGrid is only true when doing a grid of calculations */
3185  if( strcmp(chTime,"LAST") == 0 &&
3186  !(iterations.lgLastIt && !grid.lgGrid ) &&
3187  punch.lgHashEndIter[ipPun] &&
3188  punch.lg_separate_iterations[ipPun] &&
3189  !punch.lgFITS[ipPun] )
3190  {
3191  if( dynamics.lgStatic && strcmp( punch.chHashString , "TIME_DEP" )==0 )
3192  {
3193  fprintf( punch.ipPnunit[ipPun], "\"time=%f\n",
3195  }
3196  else
3197  {
3198  fprintf( punch.ipPnunit[ipPun], "%s\n",
3199  punch.chHashString );
3200  }
3201  }
3202  if( punch.lgFLUSH )
3203  fflush( punch.ipPnunit[ipPun] );
3204  }
3205  }
3206  return;
3207 }
3208 
3209 /*PunLineIntensity produce the 'punch lines intensity' output */
3210 STATIC void PunLineIntensity(FILE * ioPUN)
3211 {
3212  char chCard[INPUT_LINE_LENGTH];
3213  bool lgEOF;
3214  long int i;
3215 
3216  DEBUG_ENTRY( "PunLineIntensity()" );
3217 
3218  /* used to punch out all the emission line intensities
3219  * first initialize the line image reader */
3220 
3221  /* start the file with the input commands */
3222  input_init();
3223  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3224 
3225  lgEOF = false;
3226  while( !lgEOF )
3227  {
3228  input_readarray(chCard,&lgEOF);
3229  if( !lgEOF )
3230  {
3231  fprintf( ioPUN, "%s\n", chCard );
3232  }
3233  }
3234 
3235  /* now print any cautions or warnings */
3236  cdWarnings( ioPUN);
3237  cdCautions( ioPUN);
3238  fprintf( ioPUN, "zone=%5ld\n", nzone );
3239  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3240  fprintf( ioPUN, "begin emission lines\n" );
3241 
3242  /* >>chng 96 jul 03, only punch non-zero intensities */
3243  PunResults1Line(ioPUN," ",0,0.,"Start");
3244  for( i=0; i < LineSave.nsum; i++ )
3245  {
3246  if( LineSv[i].sumlin[LineSave.lgLineEmergent] > 0. )
3247  {
3248  PunResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
3249  LineSv[i].sumlin[LineSave.lgLineEmergent], "Line ");
3250  }
3251  }
3252 
3253  PunResults1Line(ioPUN," ",0,0.,"Flush");
3254 
3255  fprintf( ioPUN, " \n" );
3256  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3257 
3258  return;
3259 }
3260 
3261 /* lgPunchOpticalDepths true says punch optical depths */
3263 
3264 /*PunchLineStuff punch optical depths or source functions for all transferred lines */
3266  FILE * ioPUN,
3267  const char *chJob ,
3268  realnum xLimit )
3269 {
3270 
3271  long int nelem , ipLo , ipHi , i , ipISO;
3272  long index = 0;
3273  static bool lgFirst=true;
3274 
3275  DEBUG_ENTRY( "PunchLineStuff()" );
3276 
3277  /*find out which job this is and set a flag to use later */
3278  if( strcmp( &*chJob , "optical" ) == 0 )
3279  {
3280  /* punch line optical depths */
3281  lgPunchOpticalDepths = true;
3282  lgPopsFirstCall = false;
3283  }
3284  else if( strcmp( &*chJob , "populat" ) == 0 )
3285  {
3286  lgPunchOpticalDepths = false;
3287  /* level population information */
3288  if( lgFirst )
3289  {
3290  lgPopsFirstCall = true;
3291  fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
3292  lgFirst = false;
3293  }
3294  else
3295  {
3296  lgPopsFirstCall = false;
3297  }
3298  }
3299  else
3300  {
3301  fprintf( ioQQQ, " insane job in PunchLineStuff =%s\n",
3302  &*chJob );
3303  cdEXIT(EXIT_FAILURE);
3304  }
3305 
3306  /* loop over all lines, calling put1Line to create info (routine located below) */
3307  /* hydrogen like lines */
3308  /* >>chng 02 may 16, had been explicit H and He-like loops */
3309  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
3310  {
3311  for( nelem=ipISO; nelem < LIMELM; nelem++ )
3312  {
3313  if( dense.lgElmtOn[nelem] )
3314  {
3315  /* 06 aug 28, from numLevels_max to _local. */
3316  for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
3317  {
3318  for( ipLo=0; ipLo <ipHi; ipLo++ )
3319  {
3320  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
3321  continue;
3322 
3323  ++index;
3324  pun1Line( &Transitions[ipISO][nelem][ipHi][ipLo] , ioPUN , xLimit , index , dense.xIonDense[nelem][nelem+1-ipISO] );
3325  }
3326  }
3327  /* also do extra Lyman lines if optical depths are to be done,
3328  * these are line that are included only for absorption, not in the
3329  * model atoms */
3330  if( lgPunchOpticalDepths )
3331  {
3332  /* >>chng 02 aug 23, for he-like, had starting on much too high a level since
3333  * index was number of levels - caught by Adrian Turner */
3334  /* now output extra line lines, starting one above those already done above */
3335  /*for( ipHi=iso.numLevels_max[ipISO][nelem]; ipHi < iso.nLyman[ipISO]; ipHi++ )*/
3336  /* 06 aug 28, from numLevels_max to _local. */
3337  for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1].n+1; ipHi < iso.nLyman[ipISO]; ipHi++ )
3338  {
3339  ++index;
3340  pun1Line( &ExtraLymanLines[ipISO][nelem][ipHi] , ioPUN , xLimit , index , dense.xIonDense[nelem][nelem+1-ipISO]);
3341  }
3342  }
3343  }
3344  }
3345  }
3346 
3347  /* index from 1 to skip over dummy line */
3348  for( i=1; i < nLevel1; i++ )
3349  {
3350  ++index;
3351  pun1Line( &TauLines[i] , ioPUN , xLimit , index , 1. );
3352  }
3353 
3354  for( i=0; i < nWindLine; i++ )
3355  {
3356  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
3357  {
3358  ++index;
3359  pun1Line( &TauLine2[i] , ioPUN , xLimit , index , 1. );
3360  }
3361  }
3362 
3363  for( i=0; i < nUTA; i++ )
3364  {
3365  ++index;
3366  pun1Line( &UTALines[i] , ioPUN , xLimit , index , 1. );
3367  }
3368 
3369 
3370  /* do optical depths of FeII lines, if large atom is turned on */
3371  FeIIPunchLineStuff( ioPUN , xLimit , index);
3372 
3373  /* punch optical depths of H2 lines */
3374  H2_PunchLineStuff( ioPUN , xLimit , index);
3375 
3376  /* C12O16 */
3377  for( i=0; i < nCORotate; i++ )
3378  {
3379  ++index;
3380  pun1Line( &C12O16Rotate[i] , ioPUN , xLimit , index , 1. );
3381  }
3382 
3383  /* C13O16 */
3384  for( i=0; i < nCORotate; i++ )
3385  {
3386  ++index;
3387  pun1Line( &C13O16Rotate[i] , ioPUN , xLimit , index , 1. );
3388  }
3389  /*fprintf(ioPUN, "##################################\n"); */
3390  fprintf( ioPUN , "%s\n",punch.chHashString );
3391  return;
3392 }
3393 
3394 /*pun1Line called by PunchLineStuff to produce output for one line */
3395 void pun1Line( transition * t , FILE * ioPUN , realnum xLimit , long index , double abundance)
3396 {
3397 
3398  if( lgPunchOpticalDepths )
3399  {
3400  /* optical depths, no special first time, only print them */
3401  if( t->Emis->TauIn >= xLimit )
3402  {
3403  /* label like "C 4" or "H 1" */
3404  fprintf( ioPUN , "%2.2s%2.2s\t",
3406  elementnames.chIonStage[t->Hi->IonStg-1] );
3407 
3408  /* print wavelengths, either line in main printout labels,
3409  * or in various units in exponential notation - prt_wl is in prt.c */
3410  if( strcmp( punch.chConPunEnr[punch.ipConPun], "labl" )== 0 )
3411  {
3412  prt_wl( ioPUN , t->WLAng );
3413  }
3414  else
3415  {
3416  /* this converts energy in Rydbergs into any of the other units */
3417  fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t->EnergyWN * WAVNRYD)) );
3418  }
3419  /* print the optical depth */
3420  fprintf( ioPUN , "\t%.3f", t->Emis->TauIn );
3421  /* damping constant */
3422  fprintf(ioPUN, "\t%.3e",
3423  t->Emis->dampXvel / DoppVel.doppler[ t->Hi->nelem -1 ] );
3424  fprintf(ioPUN, "\n");
3425  }
3426  }
3427  else if( lgPopsFirstCall )
3428  {
3429  char chLbl[11];
3430  /* first call to line populations, print atomic parameters and indices */
3431  strcpy( chLbl, chLineLbl(t) );
3432  fprintf(ioPUN, "%li\t%s" , index , chLbl );
3433  /* stat weights */
3434  fprintf(ioPUN, "\t%.0f\t%.0f",
3435  t->Lo->g ,t->Hi->g);
3436  /* energy difference, gf */
3437  fprintf(ioPUN, "\t%.2f\t%.3e",
3438  t->EnergyWN ,t->Emis->gf);
3439  fprintf(ioPUN, "\n");
3440  }
3441  else
3442  {
3443  /* not first call, so do level populations and indices defined above */
3444  if( t->Hi->Pop > xLimit )
3445  {
3446  fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index,
3447  /* >>chng 05 may 08, add abundances, which for iso-seq species is
3448  * the density of the parent ion, for other lines, is unity.
3449  * had not been included so pops for iso seq were rel to parent ion.
3450  * caught by John Evrett */
3451  t->Lo->Pop*abundance ,
3452  t->Hi->Pop*abundance );
3453  }
3454  }
3455 }
3456 
3457 /*PunchNewContinuum produce the 'punch new continuum' output */
3458 STATIC void PunchNewContinuum(FILE * ioPUN, long ipCon )
3459 {
3460  long int ipLo, ipHi,
3461  j ,
3462  ncells;
3463 
3464  double
3465  wllo=3500. ,
3466  wlhi=7000. ,
3467  resolution = 2.;
3468 
3469  double *energy,
3470  *cont_incid,
3471  *cont_atten,
3472  *diffuse_in,
3473  *diffuse_out,
3474  *emis_lines_out,
3475  *emis_lines_in;
3476 
3477  /* get the low limit, cp_range_min is zero if not set */
3478  wllo = MAX2( rfield.anu[0] , punch.cp_range_min[ipCon] );
3479  if( punch.cp_range_max[ipCon] > 0. )
3480  {
3481  /* set to smaller of these two */
3482  wlhi = MIN2( rfield.anu[rfield.nflux-1] , punch.cp_range_max[ipCon] );
3483  }
3484  else
3485  {
3486  /* use high-energy limit since not set */
3487  wlhi = rfield.anu[rfield.nflux-1];
3488  }
3489 
3490  if( punch.cp_resolving_power[ipCon] != 0. )
3491  {
3492  /* next two bogus to fool lint - they cannot happen */
3493  ipLo = LONG_MAX;
3494  ipHi = LONG_MAX;
3495  /* we will set a new continuum mesh */
3496  ncells = (long)((wlhi-wllo)/resolution + 1);
3497  }
3498  else
3499  {
3500  /* use native continuum mesh */
3501  ipLo = ipoint(wllo)-1;
3502  ipHi = ipoint(wlhi)-1;
3503  ncells = ipHi - ipLo + 1;
3504  }
3505 
3506  /* now allocate the space */
3507  energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
3508  cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
3509  cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
3510  diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
3511  diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
3512  emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
3513  emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
3514  /*emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));
3515  emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));*/
3516 
3517  /* was the resolution changed from the default native resolution ? */
3518  if( punch.cp_resolving_power[ipCon] != 0. )
3519  {
3520  /* now create energy grid */
3521  energy[0] = wlhi;
3522  j = 0;
3523  while( energy[j]-resolution > wllo )
3524  {
3525  ++j;
3526  ASSERT( j< ncells+1 );
3527  energy[j] = energy[j-1] - resolution;
3528  }
3529 
3530  ncells = j;
3531  /* now convert to rydbergs */
3532  for( j=0; j<ncells; ++j )
3533  {
3534  energy[j] = RYDLAM / energy[j];
3535  }
3536  }
3537  else
3538  {
3539  /* now convert to rydbergs */
3540  for( j=0; j<ncells; ++j )
3541  {
3542  energy[j] = rfield.AnuOrg[j+ipLo] - rfield.widflx[j+ipLo]/2.;
3543  }
3544  }
3545 
3546  /* for cdSPEC the energy vector is the lower edge of the energy cell */
3547  /* get incident continuum */
3548  cdSPEC( 1 , energy , ncells , cont_incid );
3549  /* get attenuated incident continuum */
3550  cdSPEC( 2 , energy , ncells , cont_atten );
3551  /* get attenuated incident continuum */
3552  cdSPEC( 5 , energy , ncells , diffuse_in );
3553  /* get attenuated incident continuum */
3554  cdSPEC( 4 , energy , ncells , diffuse_out );
3556  /* different parts of lines */
3557  cdSPEC( 6 , energy , ncells , emis_lines_out );
3558  cdSPEC( 7 , energy , ncells , emis_lines_in );
3559  /*cdSPEC( 8 , energy , ncells , emis_lines_pump_out );
3560  cdSPEC( 9 , energy , ncells , emis_lines_pump_in );*/
3561 
3562  /* for this example we will do a wavelength range */
3563  for( j=0; j<ncells-1; ++j )
3564  {
3565  /* photon energy in appropriate energy or wavelength units */
3566  fprintf( ioPUN,"%.3e\t", AnuUnit((realnum)(energy[j]+rfield.widflx[j+ipLo]/2.) ) );
3567  fprintf( ioPUN,"%.3e\t", cont_incid[j] );
3568  fprintf( ioPUN,"%.3e\t", cont_atten[j] );
3569  fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
3570  fprintf( ioPUN,"%.3e",
3571  emis_lines_out[j]+emis_lines_in[j]/*+emis_lines_pump_out[j]+emis_lines_pump_in[j]*/ );
3572  fprintf( ioPUN,"\n" );
3573  }
3574 
3575  free(energy);
3576  free(cont_incid);
3577  free(diffuse_in);
3578  free(diffuse_out);
3579  free(cont_atten);
3580  free(emis_lines_out);
3581  free(emis_lines_in);
3582  /*free(emis_lines_pump_out);
3583  free(emis_lines_pump_in );*/
3584 }
3585 
3586 /* punch AGN3 hemiss, for Chapter 4, routine is below */
3587 STATIC void AGN_Hemis(FILE *ioPUN )
3588 {
3589  const int NTE = 4;
3590  realnum te[NTE] = {5000., 10000., 15000., 20000.};
3591  realnum *agn_continuum[NTE];
3592  double TempSave = phycon.te;
3593  long i , j;
3594 
3595  DEBUG_ENTRY( "AGN_Hemis()" );
3596 
3597  /* make table of continuous emission at various temperatuers */
3598  /* first allocate space */
3599  for( i=0;i<NTE; ++i)
3600  {
3601  agn_continuum[i] = (realnum *)MALLOC((unsigned)rfield.nflux*sizeof(realnum) );
3602 
3603  /* set the next temperature */
3604  /* recompute everything at this new temp */
3605  TempChange(te[i] , true);
3606  /* converge the pressure-temperature-ionization solution for this zone */
3608 
3609  /* now get the thermal emission */
3610  RT_diffuse();
3611  for(j=0;j<rfield.nflux; ++j )
3612  {
3613  agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(realnum)dense.eden/
3615  }
3616  }
3617 
3618  /* print title for line */
3619  fprintf(ioPUN,"wl");
3620  for( i=0;i<NTE; ++i)
3621  {
3622  fprintf(ioPUN,"\tT=%.0f",te[i]);
3623  }
3624  fprintf( ioPUN , "\tcont\n");
3625 
3626  /* not print all n temperatures across a line */
3627  for(j=0;j<rfield.nflux; ++j )
3628  {
3629  fprintf( ioPUN , "%12.5e",
3630  AnuUnit(rfield.AnuOrg[j]) );
3631  /* loop over the temperatures, and for each, calculate a continuum */
3632  for( i=0;i<NTE; ++i)
3633  {
3634  fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
3635  }
3636  /* cont label and end of line*/
3637  fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j]);
3638  }
3639 
3640  /* now free the continua */
3641  for( i=0;i<NTE; ++i)
3642  {
3643  free( agn_continuum[i] );
3644  }
3645 
3646  /* Restore temperature stored before this routine was called */
3647  /* and force update */
3648  TempChange(TempSave , true);
3649 
3650  fprintf( ioQQQ, "AGN_Hemis - result of punch AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
3651  cdEXIT(EXIT_FAILURE);
3652 }
3653 
3654 /*punResults punch results from punch results command */
3655 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */
3656 STATIC void punResults(FILE* ioPUN)
3657 {
3658  char chCard[INPUT_LINE_LENGTH];
3659  bool lgEOF;
3660  long int i , nelem , ion;
3661 
3662  DEBUG_ENTRY( "punResults()" );
3663 
3664  /* used to punch out line intensities, optical depths,
3665  * and column densities */
3666  lgEOF = false;
3667  /* first initialize the line image reader */
3668  input_init();
3669 
3670  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3671  lgEOF = false;
3672  input_readarray(chCard,&lgEOF);
3673  while( !lgEOF )
3674  {
3675  /* will get copy of input line image as all capital letters here */
3676  char chCAPS[INPUT_LINE_LENGTH];
3677  strcpy( chCAPS , chCard );
3678  caps( chCAPS );
3679  /* keyword HIDE means to hide the command - do not print it */
3680  if( !nMatch( "HIDE" , chCAPS ) )
3681  fprintf( ioPUN, "%s\n", chCard );
3682  input_readarray(chCard,&lgEOF);
3683  }
3684 
3685  /* first print any cautions or warnings */
3686  cdWarnings(ioPUN);
3687  cdCautions(ioPUN);
3688  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3689 
3690  fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );
3691 
3692  fprintf( ioPUN, "BEGIN EMISSION LINES\n" );
3693  PunResults1Line(ioPUN," ",0,0.,"Start");
3694 
3695  for( i=0; i < LineSave.nsum; i++ )
3696  {
3697  if( LineSv[i].sumlin[LineSave.lgLineEmergent] > 0. )
3698  {
3699  PunResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
3700  LineSv[i].sumlin[LineSave.lgLineEmergent],"Line ");
3701  }
3702  }
3703 
3704  PunResults1Line(ioPUN," ",0,0.,"Flush");
3705 
3706  fprintf( ioPUN, " \n" );
3707 
3708  fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" );
3709 
3710  /* this dumps out the whole array,*/
3711  /* following loop relies on LIMELM being 30, assert it here in case
3712  * this is ever changed */
3713  ASSERT( LIMELM == 30 );
3714  /* this order of indices is to keep 30 as the fastest variable,
3715  * and the 32 (LIMELM+1) as the slower one */
3716  for( nelem=0; nelem<LIMELM; nelem++ )
3717  {
3718  for(ion=0; ion < nelem+1; ion++)
3719  {
3720  fprintf( ioPUN, " %10.3e", mean.xIonMeans[0][nelem][ion] );
3721  /* throw line feed every 10 numbers */
3722  if( nelem==9|| nelem==19 || nelem==29 )
3723  {
3724  fprintf( ioPUN, "\n" );
3725  }
3726  }
3727  }
3728 
3729  fprintf( ioPUN, "END OF RESULTS\n" );
3730  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3731  return;
3732 }
3733 
3734 static const int LINEWIDTH = 6;
3735 
3736 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */
3737 /* the number of emission lines across one line of printout */
3739  FILE * ioPUN,
3740  /* 4 char + null string */
3741  const char *chLab,
3742  realnum wl,
3743  double xInten,
3744  /* null terminated string saying what to do */
3745  const char *chFunction)
3746 {
3747 
3748  long int i;
3749  static realnum wavelength[LINEWIDTH];
3750  static long ipLine;
3751  static double xIntenSave[LINEWIDTH];
3752  static char chLabSave[LINEWIDTH][5];
3753 
3754  DEBUG_ENTRY( "PunResults1Line()" );
3755 
3756  /* if LineWidth is changed then change format in write too */
3757 
3758  if( strcmp(chFunction,"Start") == 0 )
3759  {
3760  ipLine = 0;
3761  }
3762  else if( strcmp(chFunction,"Line ") == 0 )
3763  {
3764  /* save results in array so that they can be printed when done */
3765  wavelength[ipLine] = wl;
3766  strcpy( chLabSave[ipLine], chLab );
3767  xIntenSave[ipLine] = xInten;
3768 
3769  /* now increment the counter and then check if we have filled the line,
3770  * and so should print it */
3771  ++ipLine;
3772  /* do print now if we are in column mode (one line per line) or if we have filled up
3773  * the line */
3774  if( ( strcmp(punch.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
3775  {
3776  /* "array " is usual array 6 wide, "column" is one line per line */
3777  for( i=0; i < ipLine; i++ )
3778  {
3779  fprintf( ioPUN, " %4.4s ", chLabSave[i] );
3780  prt_wl( ioPUN, wavelength[i] );
3781  fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
3782  /* >>chng 02 apr 24, do not print type */
3783  /* single column for input into data base */
3784  if( strcmp(punch.chPunRltType,"column") == 0 )
3785  fprintf( ioPUN, "\n" );
3786  }
3787  /* only put cr if we did not just put one */
3788  if( strcmp(punch.chPunRltType,"array ") == 0 )
3789  fprintf( ioPUN, " \n" );
3790  ipLine = 0;
3791  }
3792  }
3793  else if( strcmp(chFunction,"Flush") == 0 )
3794  {
3795  if( ipLine > 0 )
3796  {
3797  /* this is an option to print many emission lines across an output line,
3798  * the array option, or a single column of numbers, the column option
3799  * that appears on the "punch results" and "punch intensity" commands
3800  */
3801  /* usual array 6 wide */
3802  for( i=0; i < ipLine; i++ )
3803  {
3804  fprintf( ioPUN, " %4.4s", chLabSave[i] );
3805  prt_wl( ioPUN, wavelength[i] );
3806  fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
3807  /* >>chng 02 apr 24, do not print type */
3808  /* single column for input into data base */
3809  if( strcmp(punch.chPunRltType,"column") == 0 )
3810  fprintf( ioPUN, "\n" );
3811  }
3812  if( strcmp(punch.chPunRltType,"array ") == 0 )
3813  fprintf( ioPUN, " \n" );
3814  }
3815  }
3816  else
3817  {
3818  fprintf( ioQQQ, " PunResults1Line called with insane request=%5.5s\n",
3819  chFunction );
3820  cdEXIT(EXIT_FAILURE);
3821  }
3822  return;
3823 }
3824 
3825 static const int NENR_GAUNT = 37;
3826 static const int NTE_GAUNT = 21;
3827 
3828 /*PunchGaunts called by punch gaunts command to output gaunt factors */
3829 STATIC void PunchGaunts(FILE* ioPUN)
3830 {
3831  long int i,
3832  charge,
3833  ite,
3834  j;
3835 
3836  realnum ener[NENR_GAUNT],
3837  ste[NTE_GAUNT],
3838  z;
3839  double tempsave;
3840  double g[NENR_GAUNT][NTE_GAUNT], gfac;
3841 
3842 
3843  DEBUG_ENTRY( "PunchGaunts()" );
3844 
3845  /* this routine is called from the PUNCH GAUNTS command
3846  * it drives the gaunt factor routine to punch gaunts over full range */
3847  tempsave = phycon.te;
3848 
3849  for( i=0; i < NTE_GAUNT; i++ )
3850  {
3851  ste[i] = 0.5f*i;
3852  }
3853 
3854  for( i=0; i < NENR_GAUNT; i++ )
3855  {
3856  ener[i] = 0.5f*i - 8.f;
3857  }
3858 
3859  for( charge = 1; charge<LIMELM; charge++ )
3860  {
3861  /* nuclear charge */
3862  z = (realnum)log10((double)charge);
3863 
3864  /* energy is log of energy */
3865  for( ite=0; ite < NTE_GAUNT; ite++ )
3866  {
3867  phycon.alogte = ste[ite];
3868  phycon.te = pow(10.,phycon.alogte);
3869 
3870  for( j=0; j < NENR_GAUNT; j++ )
3871  {
3872  gfac = cont_gaunt_calc( phycon.te, (double)charge, pow(10.,(double)ener[j]) );
3873  /*fprintf(ioQQQ, "z %.2e ener %.2e temp %.2e gfac %.3e \n",
3874  pow(10.,z), pow(10.,ener[j]), (double)phycon.te, gfac );*/
3875  g[j][ite] = gfac;
3876  }
3877  }
3878 
3879  /* now punch out the results */
3880  fprintf( ioPUN, " " );
3881  for( i=1; i <= NTE_GAUNT; i++ )
3882  {
3883  fprintf( ioPUN, "\t%6.1f", ste[i-1] );
3884  }
3885  fprintf( ioPUN, "\n" );
3886 
3887  for( j=0; j < NENR_GAUNT; j++ )
3888  {
3889  fprintf( ioPUN, "\t%6.2f", ener[j] );
3890  for( ite=0; ite < NTE_GAUNT; ite++ )
3891  {
3892  fprintf( ioPUN, "\t%6.2f", log10(g[j][ite]) );
3893  }
3894  fprintf( ioPUN, "\n" );
3895  }
3896 
3897  fprintf( ioPUN, " " );
3898  for( i=0; i < NTE_GAUNT; i++ )
3899  {
3900  fprintf( ioPUN, "\t%6.1f", ste[i] );
3901  }
3902  fprintf( ioPUN, "\n" );
3903 
3904  /* now do the same thing as triplets */
3905  fprintf( ioPUN, " " );
3906  for( i=0; i < NTE_GAUNT; i++ )
3907  {
3908  fprintf( ioPUN, "\t%6.1f", ste[i] );
3909  }
3910  fprintf( ioPUN, "\n" );
3911 
3912  for( i=0; i < NTE_GAUNT; i++ )
3913  {
3914  for( j=0; j < NENR_GAUNT; j++ )
3915  {
3916  fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", ste[i], ener[j],
3917  log10(g[j][i]) );
3918  }
3919  }
3920 
3921  fprintf( ioPUN, "Below is log(u), log(gamma**2), gff\n" );
3922  /* do the same thing as above, but this time print log(u) and log(gamma2) instead of temp and energy. */
3923  for( i=0; i < NTE_GAUNT; i++ )
3924  {
3925  for( j=0; j < NENR_GAUNT; j++ )
3926  {
3927  fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", 2.*z + log10(TE1RYD) - ste[i] , log10(TE1RYD)+ener[j]-ste[i],
3928  log10(g[j][i]) );
3929  }
3930  }
3931  fprintf( ioPUN, "end of charge = %li\n", charge);
3932  fprintf( ioPUN, "****************************\n");
3933  }
3934 
3935  phycon.te = tempsave;
3936  phycon.alogte = log10( phycon.te );
3937  return;
3938 }

Generated for cloudy by doxygen 1.8.1.2