00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include "cddefines.h"
00018 #include "cddrive.h"
00019 #include "physconst.h"
00020 #include "mean.h"
00021 #include "taulines.h"
00022 #include "struc.h"
00023 #include "iso.h"
00024 #include "mole.h"
00025 #include "hyperfine.h"
00026 #include "rt.h"
00027 #include "lines_service.h"
00028 #include "doppvel.h"
00029 #include "dense.h"
00030 #include "h2.h"
00031 #include "magnetic.h"
00032 #include "hydrogenic.h"
00033 #include "secondaries.h"
00034 #include "grainvar.h"
00035 #include "lines.h"
00036 #include "dynamics.h"
00037 #include "colden.h"
00038 #include "continuum.h"
00039 #include "ionbal.h"
00040 #include "yield.h"
00041 #include "prt.h"
00042 #include "iterations.h"
00043 #include "heavy.h"
00044 #include "conv.h"
00045 #include "geometry.h"
00046 #include "called.h"
00047 #include "helike.h"
00048 #include "opacity.h"
00049 #include "rfield.h"
00050 #include "phycon.h"
00051 #include "timesc.h"
00052 #include "radius.h"
00053 #include "atomfeii.h"
00054 #include "assertresults.h"
00055 #include "thermal.h"
00056 #include "wind.h"
00057 #include "hmi.h"
00058 #include "pressure.h"
00059 #include "elementnames.h"
00060 #include "ipoint.h"
00061 #include "gammas.h"
00062 #include "atmdat.h"
00063 #include "map.h"
00064 #include "input.h"
00065 #include "punch.h"
00066 #include "optimize.h"
00067 #include "grid.h"
00068
00069
00070
00071 static void PunResults1Line(
00072 FILE * ioPUN,
00073
00074 const char *chLab,
00075 float wl,
00076 double xInten,
00077 const char *chFunction);
00078
00079
00080 static void PunchGaunts(FILE* ioPUN);
00081
00082
00083
00084 static void punResults(FILE* ioPUN);
00085
00086 static void PunchLineStuff(
00087 FILE * ioPUN,
00088 const char *chJob ,
00089 float xLimit);
00090
00091
00092 static void AGN_Hemis(FILE *ioPUN );
00093
00094
00095 static void PunchNewContinuum(FILE * ioPUN , long ipCon );
00096
00097
00098 static void PunLineIntensity(FILE * ioPUN);
00099
00100 char *chDummy;
00101
00102 void PunchDo(
00103
00104 const char *chTime)
00105 {
00106 bool lgOK,
00107 lgTlkSav;
00108 long int
00109 ipPun,
00110 i,
00111 i1,
00112 ion,
00113 ipConMax,
00114 ipHi,
00115 ipLinMax,
00116 ipLo,
00117 ips,
00118 j,
00119 jj,
00120 limit,
00121 nd,
00122 nelem;
00123 double ConMax,
00124 RateInter,
00125 av,
00126 conem,
00127 eps,
00128 flxatt,
00129 flxcor,
00130 flxin,
00131 flxref,
00132 flxtrn,
00133 fout,
00134 fref,
00135 fsum,
00136 opConSum,
00137 opLinSum,
00138 stage,
00139 sum,
00140 texc,
00141 xLinMax;
00142
00143 DEBUG_ENTRY( "PunchDo()" );
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 if( punch.npunch < 1 )
00168 {
00169 DEBUG_EXIT( "PunchDo()" );
00170 return;
00171 }
00172
00173 for( ipPun=0; ipPun < punch.npunch; ipPun++ )
00174 {
00175
00176 punch.ipConPun = ipPun;
00177
00178
00179
00180
00181
00182
00183 if( iterations.lgLastIt || (!punch.lgPunLstIter[ipPun]) )
00184 {
00185
00186 if( strcmp(punch.chPunch[ipPun],"ABUN") == 0 )
00187 {
00188
00189 if( strcmp(chTime,"LAST") != 0 )
00190 {
00191 fprintf( punch.ipPnunit[ipPun], "%.2f",
00192 log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) );
00193 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00194 {
00195
00196
00197 fprintf( punch.ipPnunit[ipPun], "\t%.2f",
00198 log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) );
00199 }
00200 fprintf( punch.ipPnunit[ipPun], "\n" );
00201 }
00202 }
00203
00204 else if( strcmp(punch.chPunch[ipPun],"21CM") == 0 )
00205 {
00206
00207 if( strcmp(chTime,"LAST") != 0 )
00208 {
00209 fprintf( punch.ipPnunit[ipPun],
00210 "%.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",
00211
00212 radius.depth_mid_zone,
00213 hyperfine.Tspin21cm ,
00214 phycon.te ,
00215
00216 3.84e-7* EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon /
00217 SDIV( HFLines[0].TauCon ),
00218
00219 HFLines[0].PopLo ,
00220 HFLines[0].PopHi ,
00221 OccupationNumberLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ),
00222 HFLines[0].TauCon ,
00223 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon,
00224 HFLines[0].PopOpc,
00225
00226 (dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s])/SDIV( hyperfine.Tspin21cm),
00227
00228
00229 HFLines[0].TauIn,
00230 TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) ,
00231 colden.H0_ov_Tspin,
00232
00233 colden.H0_21cm_lower,
00234 colden.H0_21cm_upper,
00235 -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower)
00236 );
00237 }
00238 }
00239
00240 else if( strcmp(punch.chPunch[ipPun],"AGES") == 0 )
00241 {
00242
00243 if( strcmp(chTime,"LAST") != 0 )
00244 {
00245 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00246
00247 radius.depth_mid_zone,
00248
00249 dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot,
00250
00251 timesc.time_H2_Dest_here,
00252
00253 timesc.AgeCOMoleDest[findspecies("CO")->index],
00254
00255 timesc.AgeCOMoleDest[findspecies("OH")->index],
00256
00257 1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) );
00258 }
00259 }
00260
00261 else if( strcmp(punch.chPunch[ipPun]," AGN") == 0 )
00262 {
00263 if( strcmp(chTime,"LAST") == 0 )
00264 {
00265 if( strcmp( punch.chPunchArgs[ipPun], "HECS" ) == 0 )
00266 {
00267
00268 AGN_He1_CS(punch.ipPnunit[ipPun]);
00269 }
00270 if( strcmp( punch.chPunchArgs[ipPun], "HEMI" ) == 0 )
00271 {
00272
00273 AGN_Hemis(punch.ipPnunit[ipPun]);
00274 }
00275 else
00276 {
00277 fprintf( ioQQQ, " PunchDo does not recognize flag %4.4s set for AGN punch. This is impossible.\n",
00278 punch.chPunch[ipPun] );
00279 ShowMe();
00280 puts( "[Stop in PunchDo]" );
00281 cdEXIT(EXIT_FAILURE);
00282 }
00283 }
00284 }
00285
00286 else if( strcmp(punch.chPunch[ipPun],"ASSE") == 0 )
00287 {
00288 if( strcmp(chTime,"LAST") == 0 )
00289 {
00290
00291
00292 lgCheckAsserts(punch.ipPnunit[ipPun]);
00293
00294 }
00295 }
00296
00297 else if( strcmp(punch.chPunch[ipPun],"AVER") == 0 )
00298 {
00299 if( strcmp(chTime,"LAST") == 0 )
00300 {
00301
00302 punch_average( ipPun , "PUNCH", chDummy );
00303 }
00304 }
00305
00306 else if( strncmp(punch.chPunch[ipPun],"CHA",3) == 0 )
00307 {
00308 if( strcmp(chTime,"LAST") == 0 )
00309 {
00310
00311 ChargTranPun( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] );
00312 }
00313 }
00314
00315 else if( strcmp(punch.chPunch[ipPun],"COOL") == 0 )
00316 {
00317
00318 if( strcmp(chTime,"LAST") != 0 )
00319 CoolPunch(punch.ipPnunit[ipPun]);
00320 }
00321
00322 else if( strncmp(punch.chPunch[ipPun],"DYN" , 3) == 0 )
00323 {
00324
00325 if( strcmp(chTime,"LAST") != 0 )
00326 DynaPunch(punch.ipPnunit[ipPun] ,punch.chPunch[ipPun][3] );
00327 }
00328
00329 else if( strcmp(punch.chPunch[ipPun],"ENTH") == 0 )
00330 {
00331 if( strcmp(chTime,"LAST") != 0 )
00332 fprintf( punch.ipPnunit[ipPun],
00333 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00334 radius.depth_mid_zone,
00335 phycon.EnthalpyDensity,
00336 phycon.EnergyExcitation,
00337 phycon.EnergyIonization,
00338 phycon.EnergyBinding ,
00339 0.5*POW2(wind.windv)*dense.xMassDensity ,
00340 5./2.*pressure.PresGasCurr ,
00341 magnetic.EnthalpyDensity);
00342 }
00343
00344 else if( strcmp(punch.chPunch[ipPun],"COLU") == 0 )
00345 {
00346
00347 if( strcmp(chTime,"LAST") == 0 )
00348 {
00349 PrtColumns(punch.ipPnunit[ipPun]);
00350 }
00351 }
00352
00353 else if( strcmp(punch.chPunch[ipPun],"COLS") == 0 )
00354 {
00355
00356 if( strcmp(chTime,"LAST") == 0 )
00357 {
00358 char chHeader[2];
00359 punch_colden(chHeader , punch.ipPnunit[ipPun] , "PUNS" );
00360 }
00361 }
00362
00363 else if( strcmp(punch.chPunch[ipPun],"COMP") == 0 )
00364 {
00365
00366 if( strcmp(chTime,"LAST") != 0 )
00367 {
00368 for( jj=0; jj<rfield.nflux; jj = jj + punch.ncPunchSkip)
00369 {
00370 fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e\n",
00371 rfield.anu[jj], rfield.comup[jj]/rfield.widflx[jj],
00372 rfield.comdn[jj]/rfield.widflx[jj] );
00373 }
00374 }
00375 }
00376
00377
00378 else if( strcmp(punch.chPunch[ipPun],"CON ") == 0 )
00379 {
00380
00381
00382
00383
00384 bool lgPrintThis =false;
00385 if( punch.lgPunchEveryZone[ipPun] )
00386 {
00387
00388 if( strcmp(chTime,"LAST") != 0 )
00389 {
00390
00391 if( nzone==1 )
00392 {
00393 lgPrintThis = true;
00394 }
00395 else if( nzone%punch.nPunchEveryZone[ipPun]==0 )
00396 {
00397 lgPrintThis = true;
00398 }
00399 }
00400 else
00401 {
00402
00403 if( nzone!=1 && nzone%punch.nPunchEveryZone[ipPun]!=0 )
00404 {
00405 lgPrintThis = true;
00406 }
00407 }
00408 }
00409 else
00410 {
00411
00412 if( strcmp(chTime,"LAST") == 0 )
00413 lgPrintThis = true;
00414 }
00415 ASSERT( !punch.lgPunchEveryZone[ipPun] || punch.nPunchEveryZone[ipPun]>0 );
00416 if( lgPrintThis )
00417 {
00418 if( punch.lgPunchEveryZone[ipPun] && nzone!=1)
00419 fprintf( punch.ipPnunit[ipPun], "%s\n",
00420 punch.chHashString );
00421
00422 for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
00423 {
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434 flxin = rfield.FluxSave[j]*rfield.anu2[j]*
00435 EN1RYD/rfield.widflx[j];
00436
00437
00438 flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/rfield.widflx[j] +
00439 rfield.anu[j]*punch.PunchLWidth*rfield.reflin[j])*EN1RYD;
00440
00441
00442 flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/
00443 rfield.widflx[j]*radius.r1r0sq * rfield.trans_coef_total[j];
00444
00445
00446 conem = ((rfield.ConEmitOut[j])/
00447 rfield.widflx[j]*rfield.anu2[j] + punch.PunchLWidth*
00448 rfield.outlin[j]*rfield.anu[j])*radius.r1r0sq*
00449 EN1RYD*geometry.covgeo;
00450
00451
00452 flxtrn = conem + flxatt;
00453
00454
00455 fprintf( punch.ipPnunit[ipPun],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) );
00456
00457 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxin );
00458
00459 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxatt );
00460
00461 fprintf( punch.ipPnunit[ipPun],"%.3e\t", conem );
00462
00463 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn );
00464
00465 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref );
00466
00467 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref + flxtrn );
00468 fprintf( punch.ipPnunit[ipPun], "%4.4s\t%4.4s\t",
00469
00470 rfield.chLineLabel[j] ,
00471
00472 rfield.chContLabel[j] );
00473
00474
00475 fprintf( punch.ipPnunit[ipPun], "%.2f\n", rfield.line_count[j]/rfield.widflx[j]*rfield.anu[j] );
00476 }
00477 }
00478 }
00479
00480
00481
00482 else if( strcmp(punch.chPunch[ipPun],"CONN") == 0 )
00483 {
00484 if( strcmp(chTime,"LAST") == 0 )
00485
00486 PunchNewContinuum( punch.ipPnunit[ipPun] , (long)punch.punarg[ipPun][0]);
00487 }
00488
00489 else if( strcmp(punch.chPunch[ipPun],"CONC") == 0 )
00490 {
00491
00492
00493
00494 if( strcmp(chTime,"LAST") == 0 )
00495 {
00496
00497 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
00498 {
00499 flxin = rfield.FluxSave[j]*rfield.anu2[j]*
00500 EN1RYD/rfield.widflx[j];
00501
00502 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\n",
00503 AnuUnit(rfield.AnuOrg[j]), flxin );
00504 }
00505 }
00506 }
00507
00508 else if( strcmp(punch.chPunch[ipPun],"CONG") == 0 )
00509 {
00510
00511 if( strcmp(chTime,"LAST") == 0 )
00512 {
00513
00514 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
00515 {
00516 fsum = gv.GraphiteEmission[j]*rfield.anu2[j]*
00517 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
00518
00519 fout = gv.SilicateEmission[j]*rfield.anu2[j]*
00520 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
00521
00522
00523
00524
00525 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n",
00526 AnuUnit(rfield.AnuOrg[j]) , fsum , fout , fsum+fout );
00527 }
00528 }
00529 }
00530
00531 else if( strcmp(punch.chPunch[ipPun],"CONR") == 0 )
00532 {
00533
00534
00535
00536 if( strcmp(chTime,"LAST") == 0 )
00537 {
00538 if( geometry.lgSphere )
00539 {
00540 fprintf( punch.ipPnunit[ipPun], " Reflected continuum not predicted when SPHERE is set.\n" );
00541 fprintf( ioQQQ ,
00542 "\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" );
00543 puts( "[Stop in PunchDo]" );
00544 cdEXIT(EXIT_FAILURE);
00545 }
00546
00547 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
00548 {
00549
00550 flxref = rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/
00551 rfield.widflx[j]*EN1RYD;
00552
00553 fref = rfield.anu[j]*punch.PunchLWidth*
00554 rfield.reflin[j]*EN1RYD;
00555
00556 if( rfield.FluxSave[j] > 1e-25 )
00557 {
00558 av = rfield.ConRefIncid[j]/rfield.FluxSave[j];
00559 }
00560 else
00561 {
00562 av = 0.;
00563 }
00564
00565 fprintf( punch.ipPnunit[ipPun], "%12.5e%12.4e%12.4e%12.4e%12.4e %4.4s\n",
00566 AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref,
00567 av, rfield.chContLabel[j] );
00568 }
00569 }
00570 }
00571
00572 else if( strcmp(punch.chPunch[ipPun],"CNVE") == 0 )
00573 {
00574
00575 if( strcmp(chTime,"LAST") != 0 )
00576 {
00577 fprintf( punch.ipPnunit[ipPun],
00578 "%.5e\t%li\t%.4e\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n",
00579 radius.depth_mid_zone,
00580 conv.nPres2Ioniz,
00581 pressure.PresTotlCorrect,
00582 pressure.PresTotlCurr,
00583 (pressure.PresTotlCorrect - pressure.PresTotlCurr)*100./pressure.PresTotlCorrect,
00584 dense.EdenTrue,
00585 dense.eden,
00586 (dense.EdenTrue - dense.eden)*100./dense.EdenTrue,
00587 thermal.htot,
00588 thermal.ctot,
00589 (thermal.htot - thermal.ctot)*100./thermal.htot );
00590 }
00591 }
00592
00593 else if( strcmp(punch.chPunch[ipPun],"CONB") == 0 )
00594 {
00595
00596
00597
00598 if( strcmp(chTime,"LAST") != 0 )
00599 {
00600 for( j=0; j<rfield.nupper; j = j + punch.ncPunchSkip)
00601 {
00602 fprintf( punch.ipPnunit[ipPun], "%14.5e %14.5e %14.5e\n",
00603 AnuUnit(rfield.AnuOrg[j]), rfield.anu[j], rfield.widflx[j] );
00604 }
00605 }
00606 }
00607
00608 else if( strcmp(punch.chPunch[ipPun],"COND") == 0 )
00609 {
00610
00611
00612
00613 if( strcmp(chTime,"LAST") == 0 )
00614 {
00615
00616 for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
00617 {
00618
00619 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\n",
00620 AnuUnit(rfield.AnuOrg[j]),
00621 rfield.ConEmitLocal[nzone][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
00622 }
00623 }
00624 }
00625
00626 else if( strcmp(punch.chPunch[ipPun],"CONE") == 0 )
00627 {
00628
00629
00630
00631 if( strcmp(chTime,"LAST") == 0 )
00632 {
00633
00634 for( j=0; j<rfield.nflux;j = j +punch.ncPunchSkip)
00635 {
00636
00637 flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/
00638 rfield.widflx[j] + rfield.anu[j]*punch.PunchLWidth*
00639 rfield.reflin[j])*EN1RYD;
00640
00641
00642 conem = (rfield.ConEmitOut[j])/rfield.widflx[j]*rfield.anu2[j] +
00643 punch.PunchLWidth*rfield.outlin[j]*rfield.anu[j];
00644
00645 conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;
00646
00647
00648
00649 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n",
00650 AnuUnit(rfield.AnuOrg[j]),
00651 flxref,
00652 conem,
00653 flxref + conem,
00654 rfield.chLineLabel[j],
00655 rfield.chContLabel[j]
00656 );
00657 }
00658 }
00659 }
00660
00661
00662 else if( strcmp(punch.chPunch[ipPun],"CONf") == 0 )
00663 {
00664 if( strcmp(chTime,"LAST") == 0 )
00665 {
00666 long nu_hi , nskip;
00667 if( punch.punarg[ipPun][0] > 0. )
00668 {
00669
00670 j = ipFineCont( punch.punarg[ipPun][0] );
00671 }
00672 else
00673 {
00674 j = 0;
00675 }
00676
00677 nskip = (long)punch.punarg[ipPun][2];
00678
00679
00680 if( punch.punarg[ipPun][1]> 0. )
00681 {
00682 nu_hi = ipFineCont( punch.punarg[ipPun][1]);
00683 }
00684 else
00685 {
00686 nu_hi = rfield.nfine;
00687 }
00688 do
00689 {
00690 float sum1 = rfield.fine_opt_depth[j];
00691 float xnu = rfield.fine_anu[j+nskip/2];
00692 for( jj=1; jj<nskip; ++jj )
00693 {
00694 sum1 += rfield.fine_opt_depth[j+jj];
00695 }
00696 fprintf( punch.ipPnunit[ipPun],
00697 "%.6e\t%.3e\n",
00698 AnuUnit(xnu),
00699 sexp(sum1/nskip) );
00700 j += nskip;
00701 } while( j < nu_hi );
00702 }
00703 }
00704
00705 else if( strcmp(punch.chPunch[ipPun],"CONi") == 0 )
00706 {
00707
00708
00709
00710
00711
00712 if( strcmp(chTime,"LAST") != 0 )
00713 {
00714
00715 if( punch.punarg[ipPun][0] <= 0. )
00716 {
00717 i1 = 1;
00718 }
00719 else if( punch.punarg[ipPun][0] < 100. )
00720 {
00721 i1 = ipoint(punch.punarg[ipPun][0]);
00722 }
00723 else
00724 {
00725 i1 = (long int)punch.punarg[ipPun][0];
00726 }
00727
00728 fref = 0.;
00729 fout = 0.;
00730 fsum = 0.;
00731 sum = 0.;
00732 flxin = 0.;
00733
00734 for( j=i1-1; j < rfield.nflux; j++ )
00735 {
00736 fref += rfield.flux[j]*opac.opacity_abs[j];
00737 fout += rfield.otslin[j]*opac.opacity_abs[j];
00738 fsum += rfield.otscon[j]*opac.opacity_abs[j];
00739 sum += rfield.ConInterOut[j]*opac.opacity_abs[j];
00740 flxin += (rfield.outlin[j] + rfield.outlin_noplot[j])*opac.opacity_abs[j];
00741 }
00742 fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e%10.2e%10.2e\n",
00743 fref, fout, fsum, sum, flxin );
00744 }
00745 }
00746
00747 else if( strcmp(punch.chPunch[ipPun],"CONI") == 0 )
00748 {
00749
00750
00751
00752
00753 if( (punch.punarg[ipPun][2]>0) || (strcmp(chTime,"LAST") == 0) )
00754 {
00755
00756 bool lgPrt=false;
00757 if( punch.punarg[ipPun][2]>0 )
00758 fprintf(punch.ipPnunit[ipPun],"#punch every zone %li\n", nzone);
00759
00760
00761
00762
00763 if( punch.punarg[ipPun][0] <= 0. )
00764 {
00765 i1 = 1;
00766 }
00767 else if( punch.punarg[ipPun][0] < 100. )
00768 {
00769 i1 = ipoint(punch.punarg[ipPun][0]);
00770 }
00771 else
00772 {
00773 i1 = (long int)punch.punarg[ipPun][0];
00774 }
00775
00776 sum = 0.;
00777 for( j=i1-1; j < rfield.nflux; j++ )
00778 {
00779 flxcor = rfield.flux[j] +
00780 rfield.otslin[j] +
00781 rfield.otscon[j] +
00782 rfield.ConInterOut[j] +
00783 rfield.outlin[j] + rfield.outlin_noplot[j];
00784
00785 sum += flxcor*opac.opacity_abs[j];
00786 }
00787
00788 if( sum > 0. )
00789 {
00790 sum = 1./sum;
00791 }
00792 else
00793 {
00794 sum = 1.;
00795 }
00796
00797 fsum = 0.;
00798
00799 for( j=i1-1; j<rfield.nflux; ++j)
00800 {
00801 flxcor = rfield.flux[j] +
00802 rfield.otslin[j] +
00803 rfield.otscon[j] +
00804 rfield.ConInterOut[j]+
00805 rfield.outlin[j] + rfield.outlin_noplot[j];
00806
00807 fsum += flxcor*opac.opacity_abs[j];
00808
00809
00810
00811 RateInter = flxcor*opac.opacity_abs[j]*sum;
00812
00813
00814
00815 if( (RateInter >= punch.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) )
00816 {
00817 lgPrt = true;
00818
00819 fprintf( punch.ipPnunit[ipPun],
00820 "%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",
00821 j,
00822 AnuUnit(rfield.AnuOrg[j]),
00823 flxcor,
00824 flxcor*opac.opacity_abs[j],
00825 rfield.flux[j]/flxcor,
00826 rfield.otslin[j]/flxcor,
00827 rfield.otscon[j]/flxcor,
00828 (rfield.outlin[j] + rfield.outlin_noplot[j])/flxcor,
00829 rfield.ConInterOut[j]/flxcor,
00830 RateInter,
00831 fsum*sum,
00832 rfield.chLineLabel[j],
00833 rfield.chContLabel[j] );
00834 }
00835 }
00836 if( !lgPrt )
00837 {
00838
00839 fprintf(punch.ipPnunit[ipPun],
00840 " punchdo, the PUNCH IONIZING CONTINUUM command did not find a strong point, sum fsum were %.2e %.2e\n",
00841 sum,fsum);
00842 fprintf(punch.ipPnunit[ipPun],
00843 " punchdo, the low-frequency energy was %.5e Ryd\n",rfield.anu[i1-1]);
00844 }
00845 }
00846 }
00847
00848 else if( strcmp(punch.chPunch[ipPun],"CORA") == 0 )
00849 {
00850
00851
00852
00853
00854 if( strcmp(chTime,"LAST") == 0 )
00855 {
00856
00857 for( j=0;j<rfield.nflux;j = j + punch.ncPunchSkip)
00858 {
00859 fprintf( punch.ipPnunit[ipPun],
00860 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t",
00861 AnuUnit(rfield.AnuOrg[j]),
00862 rfield.flux[j],
00863 rfield.otslin[j],
00864 rfield.otscon[j],
00865 rfield.ConRefIncid[j],
00866 rfield.ConEmitReflec[j],
00867 rfield.ConInterOut[j],
00868 rfield.outlin[j]+rfield.outlin_noplot[j],
00869 rfield.ConEmitOut[j] ,
00870 rfield.chLineLabel[j],
00871 rfield.chContLabel[j]
00872 );
00873
00874 fprintf( punch.ipPnunit[ipPun], "%li\n", rfield.line_count[j] );
00875 }
00876 }
00877 }
00878
00879 else if( strcmp(punch.chPunch[ipPun],"CONo") == 0 )
00880 {
00881
00882
00883
00884
00885 if( strcmp(chTime,"LAST") == 0 )
00886 {
00887
00888 for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip)
00889 {
00890 fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e\n",
00891 AnuUnit(rfield.AnuOrg[j]),
00892 rfield.ConEmitOut[j]+ rfield.outlin[j] + rfield.outlin_noplot[j],
00893 (rfield.flux[j] + rfield.otscon[j] + rfield.otslin[j] +
00894 rfield.ConInterOut[j])*opac.opacity_abs[j]*
00895 rfield.anu[j] );
00896 }
00897 }
00898 }
00899
00900 else if( strcmp(punch.chPunch[ipPun],"CONO") == 0 )
00901 {
00902
00903
00904
00905
00906 if( strcmp(chTime,"LAST") == 0 )
00907 {
00908
00909 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
00910 {
00911 fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e%10.2e%10.2e\n",
00912 AnuUnit(rfield.AnuOrg[j]),
00913 rfield.flux[j]*rfield.anu2[j]* EN1RYD/rfield.widflx[j]*radius.r1r0sq,
00914 rfield.ConInterOut[j]/rfield.widflx[j]*rfield.anu2[j]*radius.r1r0sq*EN1RYD,
00915 punch.PunchLWidth* (rfield.outlin[j]+rfield.outlin_noplot[j])*rfield.anu[j]*radius.r1r0sq*EN1RYD,
00916 (rfield.ConInterOut[j]/rfield.widflx[j]*
00917 rfield.anu2[j] + punch.PunchLWidth*(rfield.outlin[j]+rfield.outlin_noplot[j])*
00918 rfield.anu[j])*radius.r1r0sq*EN1RYD );
00919 }
00920 }
00921 }
00922
00923 else if( strcmp(punch.chPunch[ipPun],"CONT") == 0 )
00924 {
00925
00926
00927
00928
00929
00930 if( strcmp(chTime,"LAST") == 0 )
00931 {
00932
00933 for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip)
00934 {
00935
00936
00937
00938
00939 flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/
00940 rfield.widflx[j]*radius.r1r0sq * rfield.trans_coef_total[j];
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957 conem = (rfield.ConEmitOut[j] + rfield.outlin[j]) / rfield.widflx[j]*
00958 rfield.anu2[j]*radius.r1r0sq*EN1RYD*geometry.covgeo;
00959
00960 flxtrn = conem + flxatt;
00961
00962
00963
00964
00965 fprintf( punch.ipPnunit[ipPun],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) );
00966 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn);
00967 fprintf( punch.ipPnunit[ipPun],"%.3e\n", rfield.trans_coef_total[j] );
00968 }
00969 }
00970 }
00971
00972 else if( strcmp(punch.chPunch[ipPun],"CON2") == 0 )
00973 {
00974
00975 if( strcmp(chTime,"LAST") == 0 )
00976 {
00977
00978 for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
00979 {
00980 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\n",
00981 AnuUnit(rfield.AnuOrg[j]),
00982 rfield.TotDiff2Pht[j]/rfield.widflx[j] ,
00983 rfield.TotDiff2Pht[j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
00984 }
00985 }
00986 }
00987
00988 else if( strcmp(punch.chPunch[ipPun],"DUSE") == 0 )
00989 {
00990
00991 if( strcmp(chTime,"LAST") != 0 )
00992 {
00993 fprintf( punch.ipPnunit[ipPun], " %.5e\t",
00994 radius.depth_mid_zone );
00995
00996
00997 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
00998
00999
01000 fprintf( punch.ipPnunit[ipPun], "%.2e\n" , rfield.extin_mag_V_point);
01001 }
01002 }
01003
01004 else if( strcmp(punch.chPunch[ipPun],"DUSO") == 0 )
01005 {
01006
01007 if( strcmp(chTime,"LAST") == 0 )
01008 {
01009 for( j=0; j < rfield.nflux; j++ )
01010 {
01011 double scat;
01012 fprintf( punch.ipPnunit[ipPun],
01013 "%.5e\t%.2e\t%.2e\t%.2e\t",
01014
01015 AnuUnit(rfield.AnuOrg[j]),
01016
01017 gv.dstab[j] + gv.dstsc[j],
01018
01019 gv.dstab[j],
01020
01021 gv.dstsc[j] );
01022
01023 scat = 0.;
01024
01025 for( nd=0; nd < gv.nBin; nd++ )
01026 {
01027 scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund;
01028 }
01029
01030 fprintf( punch.ipPnunit[ipPun],
01031 "%.2e", scat );
01032 fprintf( punch.ipPnunit[ipPun],
01033 "%.2e\n", gv.dstsc[j]/(gv.dstab[j] + gv.dstsc[j]) );
01034 }
01035 }
01036 }
01037
01038 else if( strcmp(punch.chPunch[ipPun],"DUSP") == 0 )
01039 {
01040
01041 if( strcmp(chTime,"LAST") != 0 )
01042 {
01043
01044 static bool lgMustPrtHeader = true;
01045
01046 if( lgMustPrtHeader )
01047 {
01048
01049 fprintf( punch.ipPnunit[ipPun], "#Depth" );
01050 for( nd=0; nd < gv.nBin; ++nd )
01051 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01052 fprintf( punch.ipPnunit[ipPun], "\n" );
01053
01054
01055 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
01056 for( nd=0; nd < gv.nBin; ++nd )
01057 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01058 fprintf( punch.ipPnunit[ipPun], "\n" );
01059
01060
01061 lgMustPrtHeader = false;
01062 }
01063 fprintf( punch.ipPnunit[ipPun], " %.5e",
01064 radius.depth_mid_zone );
01065
01066 for( nd=0; nd < gv.nBin; ++nd )
01067 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD );
01068 fprintf( punch.ipPnunit[ipPun], "\n" );
01069 }
01070 }
01071
01072 else if( strcmp(punch.chPunch[ipPun],"DUSR") == 0 )
01073 {
01074
01075 if( strcmp(chTime,"LAST") != 0 )
01076 {
01077
01078 static bool lgMustPrtHeader = true;
01079
01080 if( lgMustPrtHeader )
01081 {
01082
01083 fprintf( punch.ipPnunit[ipPun], "#Depth" );
01084 for( nd=0; nd < gv.nBin; ++nd )
01085 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01086 fprintf( punch.ipPnunit[ipPun], "\n" );
01087
01088
01089 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
01090 for( nd=0; nd < gv.nBin; ++nd )
01091 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01092 fprintf( punch.ipPnunit[ipPun], "\n" );
01093
01094
01095 lgMustPrtHeader = false;
01096 }
01097 fprintf( punch.ipPnunit[ipPun], " %.5e",
01098 radius.depth_mid_zone );
01099
01100 for( nd=0; nd < gv.nBin; ++nd )
01101 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used );
01102 fprintf( punch.ipPnunit[ipPun], "\n" );
01103 }
01104 }
01105
01106 else if( strcmp(punch.chPunch[ipPun],"DUST") == 0 )
01107 {
01108
01109 if( strcmp(chTime,"LAST") != 0 )
01110 {
01111
01112 static bool lgMustPrtHeader = true;
01113
01114 if( lgMustPrtHeader )
01115 {
01116
01117 fprintf( punch.ipPnunit[ipPun], "#Depth" );
01118 for( nd=0; nd < gv.nBin; ++nd )
01119 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01120 fprintf( punch.ipPnunit[ipPun], "\n" );
01121
01122
01123 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
01124 for( nd=0; nd < gv.nBin; ++nd )
01125 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01126 fprintf( punch.ipPnunit[ipPun], "\n" );
01127
01128
01129 lgMustPrtHeader = false;
01130 }
01131 fprintf( punch.ipPnunit[ipPun], " %.5e",
01132 radius.depth_mid_zone );
01133 for( nd=0; nd < gv.nBin; ++nd )
01134 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust );
01135 fprintf( punch.ipPnunit[ipPun], "\n" );
01136 }
01137 }
01138
01139 else if( strcmp(punch.chPunch[ipPun],"DUSC") == 0 )
01140 {
01141
01142
01143 if( strcmp(chTime,"LAST") != 0 )
01144 {
01145
01146 static bool lgMustPrtHeader = true;
01147
01148 if( lgMustPrtHeader )
01149 {
01150
01151 fprintf( punch.ipPnunit[ipPun], "#Depth\tne(grn)" );
01152 for( nd=0; nd < gv.nBin; ++nd )
01153 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01154 fprintf( punch.ipPnunit[ipPun], "\n" );
01155
01156
01157 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)\tne\t" );
01158 for( nd=0; nd < gv.nBin; ++nd )
01159 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01160 fprintf( punch.ipPnunit[ipPun], "\n" );
01161
01162
01163 lgMustPrtHeader = false;
01164 }
01165
01166 fprintf( punch.ipPnunit[ipPun], " %.5e\t%.4e",
01167 radius.depth_mid_zone ,
01168
01169
01170 gv.TotalEden );
01171
01172
01173 for( nd=0; nd < gv.nBin; ++nd )
01174 {
01175 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ );
01176 }
01177 fprintf( punch.ipPnunit[ipPun], "\n" );
01178 }
01179 }
01180
01181 else if( strcmp(punch.chPunch[ipPun],"DUSH") == 0 )
01182 {
01183
01184 if( strcmp(chTime,"LAST") != 0 )
01185 {
01186
01187 static bool lgMustPrtHeader = true;
01188
01189 if( lgMustPrtHeader )
01190 {
01191
01192 fprintf( punch.ipPnunit[ipPun], "#Depth" );
01193 for( nd=0; nd < gv.nBin; ++nd )
01194 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01195 fprintf( punch.ipPnunit[ipPun], "\n" );
01196
01197
01198 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
01199 for( nd=0; nd < gv.nBin; ++nd )
01200 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01201 fprintf( punch.ipPnunit[ipPun], "\n" );
01202
01203
01204 lgMustPrtHeader = false;
01205 }
01206 fprintf( punch.ipPnunit[ipPun], " %.5e",
01207 radius.depth_mid_zone );
01208
01209 for( nd=0; nd < gv.nBin; ++nd )
01210 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl );
01211 fprintf( punch.ipPnunit[ipPun], "\n" );
01212 }
01213 }
01214
01215 else if( strcmp(punch.chPunch[ipPun],"DUSV") == 0 )
01216 {
01217
01218 if( strcmp(chTime,"LAST") != 0 )
01219 {
01220
01221 static bool lgMustPrtHeader = true;
01222
01223 if( lgMustPrtHeader )
01224 {
01225
01226 fprintf( punch.ipPnunit[ipPun], "#Depth" );
01227 for( nd=0; nd < gv.nBin; ++nd )
01228 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
01229 fprintf( punch.ipPnunit[ipPun], "\n" );
01230
01231
01232 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
01233 for( nd=0; nd < gv.nBin; ++nd )
01234 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
01235 fprintf( punch.ipPnunit[ipPun], "\n" );
01236
01237
01238 lgMustPrtHeader = false;
01239 }
01240 fprintf( punch.ipPnunit[ipPun], " %.5e",
01241 radius.depth_mid_zone );
01242
01243 for( nd=0; nd < gv.nBin; ++nd )
01244 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 );
01245 fprintf( punch.ipPnunit[ipPun], "\n" );
01246 }
01247 }
01248
01249
01250 else if( strcmp(punch.chPunch[ipPun],"DUSQ") == 0 )
01251 {
01252
01253 if( strcmp(chTime,"LAST") == 0 )
01254 {
01255 for( j=0; j < rfield.nflux; j++ )
01256 {
01257 fprintf( punch.ipPnunit[ipPun], "%11.4e",
01258 rfield.anu[j] );
01259 for( nd=0; nd < gv.nBin; nd++ )
01260 {
01261 fprintf( punch.ipPnunit[ipPun], "%9.1e%9.1e",
01262 gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea,
01263 gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea );
01264 }
01265 fprintf( punch.ipPnunit[ipPun], "\n" );
01266 }
01267 }
01268 }
01269
01270 else if( strcmp(punch.chPunch[ipPun],"ELEM") == 0 )
01271 {
01272 if( strcmp(chTime,"LAST") != 0 )
01273 {
01274 float renorm = 1.f;
01275
01276
01277
01278 nelem = (long int)punch.punarg[ipPun][0];
01279 ASSERT( nelem >= ipHYDROGEN );
01280
01281
01282 if( dense.lgElmtOn[nelem] )
01283 {
01284
01285
01286 if( punch.punarg[ipPun][1] == 0 )
01287 renorm = dense.gas_phase[nelem];
01288
01289 fprintf( punch.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone );
01290
01291 for( j=0; j <= (nelem + 1); ++j)
01292 {
01293 fprintf( punch.ipPnunit[ipPun], "\t%.2e",
01294 dense.xIonDense[nelem][j]/renorm );
01295 }
01296 if( nelem==ipHYDROGEN )
01297 {
01298
01299 fprintf( punch.ipPnunit[ipPun], "\t%.2e",
01300 hmi.H2_total/renorm );
01301 }
01302
01303 else if( nelem==ipCARBON )
01304 {
01305 fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
01306 colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm);
01307 fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e",
01308 colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm);
01309 fprintf( punch.ipPnunit[ipPun], "\t%.2e",
01310 findspecies("CO")->hevmol/renorm );
01311 }
01312 else if( nelem==ipOXYGEN )
01313 {
01314 fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
01315 colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm);
01316 }
01317 fprintf( punch.ipPnunit[ipPun], "\n" );
01318 }
01319 }
01320 }
01321
01322 else if( strcmp(punch.chPunch[ipPun],"RECA") == 0 )
01323 {
01324
01325
01326 ion_recombAGN( punch.ipPnunit[ipPun] );
01327 cdEXIT(EXIT_FAILURE);
01328 }
01329
01330 else if( strcmp(punch.chPunch[ipPun],"RECE") == 0 )
01331 {
01332
01333
01334
01335
01336
01337
01338 fprintf( punch.ipPnunit[ipPun],
01339 "%12.4e %12.4e %12.4e %12.4e\n",
01340 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][0][ipRecRad],
01341 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][0][ipRecNetEsc] ,
01342 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecRad],
01343 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecNetEsc]);
01344 }
01345
01346 else if( strcmp(punch.chPunch[ipPun],"FEco") == 0 )
01347 {
01348
01349 if( strcmp(chTime,"LAST") == 0 )
01350 {
01351 for( j=0; j < nFeIIConBins; j++ )
01352 {
01353
01354
01355
01356
01357
01358
01361 fprintf( punch.ipPnunit[ipPun], "%.2f\t%e \n",
01362 (FeII_Cont[j][1]+FeII_Cont[j][2])/2.,
01363 FeII_Cont[j][0] );
01364 }
01365 }
01366 }
01367
01368
01369 else if( strcmp(punch.chPunch[ipPun],"FEcl") == 0 )
01370 {
01371 if( strcmp(chTime,"LAST") == 0 )
01372 {
01373
01374 FeIIPunchColden( punch.ipPnunit[ipPun] );
01375 }
01376 }
01377
01378 else if( strcmp(punch.chPunch[ipPun],"FE2l") == 0 )
01379 {
01380 if( strcmp(chTime,"LAST") == 0 )
01381 {
01382
01383 FeIIPunchLevels( punch.ipPnunit[ipPun] );
01384 }
01385 }
01386
01387 else if( strcmp(punch.chPunch[ipPun],"FEli") == 0 )
01388 {
01389 if( strcmp(chTime,"LAST") == 0 )
01390 {
01391
01392 FeIIPunchLines( punch.ipPnunit[ipPun] );
01393 }
01394 }
01395
01396 else if( strcmp(punch.chPunch[ipPun],"FE2o") == 0 )
01397 {
01398 if( strcmp(chTime,"LAST") == 0 )
01399 {
01400
01401 FeIIPunchOpticalDepth( punch.ipPnunit[ipPun] );
01402 }
01403 }
01404
01405 else if( strcmp(punch.chPunch[ipPun],"FE2f") == 0 )
01406 {
01407
01408
01409 if( strcmp(chTime,"LAST") != 0 )
01410 pfeii(punch.ipPnunit[ipPun]);
01411 }
01412
01413 else if( strcmp(punch.chPunch[ipPun],"FE2d") == 0 )
01414 {
01415
01416 if( strcmp(chTime,"LAST") != 0 )
01417 FeIIPunDepart(punch.ipPnunit[ipPun],false);
01418 }
01419
01420 else if( strcmp(punch.chPunch[ipPun],"FE2D") == 0 )
01421 {
01422
01423 if( strcmp(chTime,"LAST") != 0 )
01424 FeIIPunDepart(punch.ipPnunit[ipPun],true);
01425 }
01426
01427 else if( strcmp(punch.chPunch[ipPun],"FE2p") == 0 )
01428 {
01429 bool lgFlag = false;
01430 if( punch.punarg[ipPun][2] )
01431 lgFlag = true;
01432
01433 if( strcmp(chTime,"LAST") != 0 )
01434 FeIIPunPop(punch.ipPnunit[ipPun],false,0,0,
01435 lgFlag);
01436 }
01437
01438 else if( strcmp(punch.chPunch[ipPun],"FE2P") == 0 )
01439 {
01440 bool lgFlag = false;
01441 if( punch.punarg[ipPun][2] )
01442 lgFlag = true;
01443
01444 if( strcmp(chTime,"LAST") != 0 )
01445 FeIIPunPop(punch.ipPnunit[ipPun],
01446 true,
01447 (long int)punch.punarg[ipPun][0],
01448 (long int)punch.punarg[ipPun][1],
01449 lgFlag);
01450 }
01451
01452
01453 else if( strcmp(punch.chPunch[ipPun],"FITS") == 0 )
01454 {
01455 if( strcmp(chTime,"LAST") == 0 )
01456 {
01457 punchFITSfile( punch.ipPnunit[ipPun], NUM_OUTPUT_TYPES );
01458 }
01459 }
01460
01461 else if( strcmp(punch.chPunch[ipPun],"GAMt") == 0 )
01462 {
01463 if( strcmp(chTime,"LAST") != 0 )
01464 {
01465 long ns;
01466
01467 for( nelem=0; nelem < LIMELM; nelem++ )
01468 {
01469 for( ion=0; ion <= nelem; ion++ )
01470 {
01471 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
01472 {
01473 fprintf( punch.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e",
01474 nelem+1, ion+1, ns+1,
01475 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
01476 ionbal.PhotoRate_Shell[nelem][ion][ns][1] ,
01477 ionbal.PhotoRate_Shell[nelem][ion][ns][2] );
01478
01479 for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ )
01480 {
01481 fprintf( punch.ipPnunit[ipPun], "%5.2f",
01482 t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) );
01483 }
01484 fprintf( punch.ipPnunit[ipPun], "\n" );
01485 }
01486 }
01487 }
01488 }
01489 }
01490
01491
01492 else if( strcmp(punch.chPunch[ipPun],"GAMe") == 0 )
01493 {
01494 if( strcmp(chTime,"LAST") != 0 )
01495 {
01496 int ns;
01497 nelem = (long)punch.punarg[ipPun][0];
01498 ion = (long)punch.punarg[ipPun][1];
01499
01500 ns = Heavy.nsShells[nelem][ion]-1;
01501
01502 GammaPrt(
01503 opac.ipElement[nelem][ion][ns][0] ,
01504 opac.ipElement[nelem][ion][ns][1] ,
01505 opac.ipElement[nelem][ion][ns][2] ,
01506 punch.ipPnunit[ipPun],
01507 ionbal.PhotoRate_Shell[nelem][ion][ns][0] ,
01508 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 );
01509 }
01510 }
01511
01512 else if( strcmp(punch.chPunch[ipPun],"GAUN") == 0 )
01513 {
01514
01515 if( strcmp(chTime,"LAST") != 0 )
01516 PunchGaunts(punch.ipPnunit[ipPun]);
01517 }
01518
01519
01520 else if( strcmp(punch.chPunch[ipPun],"GRID") == 0 )
01521 {
01522 if( (strcmp(chTime,"LAST") == 0 ) && grid.lgGridDone )
01523 {
01526
01527 fprintf(punch.ipPnunit[ipPun], "PARAM1" );
01528 for( i=2; i<= grid.nintparm; i++ )
01529 {
01530 fprintf(punch.ipPnunit[ipPun], "\tPARAM%li", i );
01531 }
01532 fprintf(punch.ipPnunit[ipPun], "\n" );
01533 for( i=0; i<grid.totNumModels; i++ )
01534 {
01535 for( j=0; j< grid.nintparm; j++ )
01536 {
01537 fprintf(punch.ipPnunit[ipPun], "%f\t", grid.interpParameters[i][j] );
01538 }
01539 fprintf(punch.ipPnunit[ipPun], "\n" );
01540 }
01541 }
01542 }
01543
01544 else if( strcmp(punch.chPunch[ipPun],"HISp") == 0 )
01545 {
01546
01547 if( strcmp(chTime,"LAST") != 0 )
01548 {
01549
01550 if( !conv.lgConvPres )
01551 {
01552 fprintf( punch.ipPnunit[ipPun],
01553 "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n",
01554 iteration , nzone );
01555 }
01556
01557 if( !conv.lgConvTemp )
01558 {
01559 fprintf( punch.ipPnunit[ipPun],
01560 "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n",
01561 iteration , nzone );
01562 }
01563 for(i=0; i<conv.hist_pres_npres; ++i )
01564 {
01565
01566 fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
01567 iteration,
01568 nzone,
01569 conv.hist_pres_density[i],
01570 conv.hist_pres_current[i],
01571 conv.hist_pres_correct[i]);
01572 }
01573 }
01574 }
01575
01576 else if( strcmp(punch.chPunch[ipPun],"HISt") == 0 )
01577 {
01578
01579 if( strcmp(chTime,"LAST") != 0 )
01580 {
01581
01582 if( !conv.lgConvPres )
01583 {
01584 fprintf( punch.ipPnunit[ipPun],
01585 "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n",
01586 iteration , nzone );
01587 }
01588
01589 if( !conv.lgConvTemp )
01590 {
01591 fprintf( punch.ipPnunit[ipPun],
01592 "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n",
01593 iteration , nzone );
01594 }
01595 for(i=0; i<conv.hist_temp_ntemp; ++i )
01596 {
01597
01598 fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
01599 iteration,
01600 nzone,
01601 conv.hist_temp_temp[i],
01602 conv.hist_temp_heat[i],
01603 conv.hist_temp_cool[i]);
01604 }
01605 }
01606 }
01607
01608 else if( strncmp(punch.chPunch[ipPun],"H2",2) == 0 )
01609 {
01610
01611 H2_PunchDo( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] , chTime, ipPun );
01612 }
01613
01614 else if( strcmp(punch.chPunch[ipPun],"HEAT") == 0 )
01615 {
01616
01617 if( strcmp(chTime,"LAST") != 0 )
01618 HeatPunch(punch.ipPnunit[ipPun]);
01619 }
01620
01621
01622 else if( strcmp(punch.chPunch[ipPun],"HUMM") == 0 )
01623 {
01624 eps = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul/
01625 (EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL *dense.eden);
01626 fprintf( punch.ipPnunit[ipPun],
01627 " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
01628 radius.depth_mid_zone,
01629 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn,
01630 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][1],
01631 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]*dense.xIonDense[ipHYDROGEN][1],
01632 phycon.te, EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].damp, eps );
01633 }
01634
01635 else if( strcmp(punch.chPunch[ipPun],"HYDc") == 0 )
01636 {
01637
01638 fprintf( punch.ipPnunit[ipPun],
01639 " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01640 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden,
01641 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
01642 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
01643 hmi.H2_total/dense.gas_phase[ipHYDROGEN],
01644 hmi.Hmolec[ipMH2p]/dense.gas_phase[ipHYDROGEN],
01645 hmi.Hmolec[ipMH3p]/dense.gas_phase[ipHYDROGEN],
01646 hmi.Hmolec[ipMHm]/dense.gas_phase[ipHYDROGEN] );
01647 }
01648
01649 else if( strcmp(punch.chPunch[ipPun],"HYDi") == 0 )
01650 {
01651 if( strcmp(chTime,"LAST") != 0 )
01652 {
01653
01654
01655 RateInter = 0.;
01656 stage = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][0]*dense.eden*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
01657 fref = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
01658 fout = iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
01659
01660 for( ion=ipH2s; ion < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ion++ )
01661 {
01662
01663 RateInter +=
01664 EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Aul*
01665 (EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Pesc +
01666 EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Pelec_esc +
01667 EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Pdest);
01668
01669 fref += iso.gamnc[ipH_LIKE][ipHYDROGEN][ion]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ion];
01670
01671 stage += iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ion]*dense.eden*
01672 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ion];
01673 fout += iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ion];
01674 }
01675 fprintf( punch.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e",
01676 nzone,
01677
01678 iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s],
01679 iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s]* dense.EdenHCorr,
01680 ionbal.RateRecomTot[ipHYDROGEN][0] );
01681
01682 fprintf( punch.ipPnunit[ipPun], "\t%.2e",
01683 iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN] );
01684
01685 fprintf( punch.ipPnunit[ipPun],
01686 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01687 dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0],
01688 iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/(ionbal.RateRecomTot[ipHYDROGEN][0]),
01689 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][1][ipRecEsc],
01690 RateInter,
01691 fref/MAX2(1e-37,fout),
01692 stage/MAX2(1e-37,fout),
01693
01694 iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][0]/(dense.eden* dense.xIonDense[ipHYDROGEN][1]) ,
01695 secondaries.csupra[ipHYDROGEN][0]);
01696
01697 GammaPrt(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0],rfield.nflux,iso.ipOpac[ipH_LIKE][ipHYDROGEN][ipH1s],
01698 punch.ipPnunit[ipPun],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*
01699 0.05);
01700 }
01701 }
01702
01703 else if( strcmp(punch.chPunch[ipPun],"HYDp") == 0 )
01704 {
01705
01706 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e",
01707 radius.depth_mid_zone,
01708 dense.xIonDense[ipHYDROGEN][0],
01709 dense.xIonDense[ipHYDROGEN][1] );
01710
01711 for( j=ipH1s; j < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]-1; j++ )
01712 {
01713 fprintf( punch.ipPnunit[ipPun], "\t%.2e",
01714 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][j]*dense.xIonDense[ipHYDROGEN][1] );
01715 }
01716 fprintf( punch.ipPnunit[ipPun], "\n" );
01717 }
01718
01719 else if( strcmp(punch.chPunch[ipPun],"HYD2") == 0 )
01720 {
01721
01722 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\n",
01723 radius.depth_mid_zone, HFLines[0].TauIn , phycon.te );
01724 }
01725
01726 else if( strcmp(punch.chPunch[ipPun],"HYDl") == 0 )
01727 {
01728 if( strcmp(chTime,"LAST") == 0 )
01729 {
01730
01731
01732 for( ipHi=4; ipHi<iso.numLevels_local[ipH_LIKE][ipHYDROGEN]-1; ++ipHi )
01733 {
01734 for( ipLo=ipHi-5; ipLo<ipHi; ++ipLo )
01735 {
01736 if( ipLo<0 )
01737 continue;
01738 fprintf(punch.ipPnunit[ipPun], "%li\t%li\t%.4e\t%.2e\n",
01739 ipHi, ipLo,
01740
01741 1./EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].EnergyWN,
01742 EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].TauIn );
01743 }
01744 }
01745 }
01746 }
01747
01748
01749 else if( strcmp(punch.chPunch[ipPun],"HYDL") == 0 )
01750 {
01751
01752 double popul = iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]/SDIV(iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]);
01753
01754 texc = TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );
01755 fprintf( punch.ipPnunit[ipPun],
01756 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
01757 radius.depth_mid_zone,
01758 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn,
01759 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot,
01760 popul,
01761 texc,
01762 phycon.te,
01763 texc/phycon.te ,
01764 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc,
01765 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest,
01766 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].pump,
01767 opac.opacity_abs[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1],
01768 opac.albedo[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] );
01769 }
01770
01771 else if( strcmp(punch.chPunch[ipPun],"HYDr") == 0 )
01772 {
01773
01774 phycon.te = 2500.;
01775 tfidle(false);
01776 while( phycon.te <= 20000. )
01777 {
01778 double r1;
01779 double ThinCoolingCaseB;
01780
01781 tfidle(false);
01782
01783 r1 = HydroRecCool(1,0);
01784 ThinCoolingCaseB = pow(10.,((-25.859117 +
01785 0.16229407*phycon.telogn[0] +
01786 0.34912863*phycon.telogn[1] -
01787 0.10615964*phycon.telogn[2])/(1. +
01788 0.050866793*phycon.telogn[0] -
01789 0.014118924*phycon.telogn[1] +
01790 0.0044980897*phycon.telogn[2] +
01791 6.0969594e-5*phycon.telogn[3])))/phycon.te;
01792
01793 fprintf( punch.ipPnunit[ipPun], " %10.2e\t",
01794 phycon.te);
01795 fprintf( punch.ipPnunit[ipPun], " %10.2e\t",
01796 (r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) );
01797
01798 fprintf( punch.ipPnunit[ipPun], " %10.2e\t",
01799 r1/(BOLTZMANN*phycon.te));
01800
01801 fprintf( punch.ipPnunit[ipPun], " %10.2e\n",
01802 ThinCoolingCaseB/(BOLTZMANN*phycon.te));
01803
01804 phycon.te = phycon.te *2.f;
01805 tfidle(false);
01806 }
01807 }
01808
01809 else if( strcmp(punch.chPunch[ipPun],"IONI") == 0 )
01810 {
01811 if( strcmp(chTime,"LAST") == 0 )
01812 {
01813
01814 PrtMeanIon( 'i', false , punch.ipPnunit[ipPun] );
01815 }
01816 }
01817
01818
01819 else if( strcmp(punch.chPunch[ipPun],"IONR") == 0 )
01820 {
01821 if( strcmp(chTime,"LAST") != 0 )
01822 {
01823
01824 nelem = (long)punch.punarg[ipPun][0];
01825 fprintf( punch.ipPnunit[ipPun],
01826 "%.5e\t%.4e\t%.4e",
01827 radius.depth_mid_zone,
01828 dense.eden ,
01829 dynamics.Rate);
01830
01831
01832 for( ion=0; ion<nelem+1; ++ion )
01833 {
01834 fprintf( punch.ipPnunit[ipPun],
01835 "\t%.4e\t%.4e\t%.4e\t%.4e",
01836 dense.xIonDense[nelem][ion] ,
01837 ionbal.RateIonizTot[nelem][ion] ,
01838 ionbal.RateRecomTot[nelem][ion] ,
01839 dynamics.Source[nelem][ion] );
01840 }
01841 fprintf( punch.ipPnunit[ipPun], "\n");
01842 }
01843 }
01844
01845 else if( strcmp(punch.chPunch[ipPun]," IP ") == 0 )
01846 {
01847 if( strcmp(chTime,"LAST") == 0 )
01848 {
01849
01850 for( nelem=0; nelem < LIMELM; nelem++ )
01851 {
01852 int ion_big;
01853 double energy;
01854
01855
01856 const int NELEM_LINE = 10;
01857
01858 for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE )
01859 {
01860 int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem);
01861
01862
01863 fprintf( punch.ipPnunit[ipPun],
01864 "\n%2.2s", elementnames.chElementSym[nelem]);
01865
01866
01867 for( ion=ion_big; ion <= ion_limit; ++ion )
01868 {
01869 fprintf( punch.ipPnunit[ipPun], "\t%4ld", ion+1 );
01870 }
01871 fprintf( punch.ipPnunit[ipPun], "\n" );
01872
01873
01874 ASSERT( ion_limit < LIMELM );
01875
01876 for( ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ )
01877 {
01878
01879
01880 fprintf( punch.ipPnunit[ipPun], "%2.2s", Heavy.chShell[ips]);
01881
01882
01883 for( ion=ion_big; ion<=ion_limit; ++ion )
01884 {
01885
01886
01887
01888 if( ips >= Heavy.nsShells[nelem][ion] )
01889 break;
01890
01891
01892 energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0);
01893
01894
01895 if( energy < 10. )
01896 {
01897 fprintf( punch.ipPnunit[ipPun], "\t%6.3f", energy );
01898 }
01899 else if( energy < 100. )
01900 {
01901 fprintf( punch.ipPnunit[ipPun], "\t%6.2f", energy );
01902 }
01903 else if( energy < 1000. )
01904 {
01905 fprintf( punch.ipPnunit[ipPun], "\t%6.1f", energy );
01906 }
01907 else
01908 {
01909 fprintf( punch.ipPnunit[ipPun], "\t%6ld", (long)(energy) );
01910 }
01911 }
01912
01913
01914 fprintf( punch.ipPnunit[ipPun], "\n" );
01915 }
01916 }
01917 }
01918 }
01919 }
01920
01921 else if( strcmp(punch.chPunch[ipPun],"LINC") == 0 )
01922 {
01923
01924
01925 if( strcmp(chTime,"LAST") != 0 )
01926 {
01927
01928 lgOK = true;
01929 punch_line(punch.ipPnunit[ipPun],"PUNC",lgOK,chDummy);
01930 }
01931 }
01932
01933 else if( strcmp(punch.chPunch[ipPun],"LIND") == 0 )
01934 {
01935
01936 PunchLineData(punch.ipPnunit[ipPun]);
01937
01938
01939
01940
01941 }
01942
01943 else if( strcmp(punch.chPunch[ipPun],"LINL") == 0 )
01944 {
01945
01946 prt_LineLabels(punch.ipPnunit[ipPun]);
01947 }
01948
01949 else if( strcmp(punch.chPunch[ipPun],"LINO") == 0 )
01950 {
01951 if( strcmp(chTime,"LAST") == 0 )
01952 {
01953
01954 PunchLineStuff(punch.ipPnunit[ipPun],"optical" , punch.punarg[ipPun][0]);
01955 }
01956 }
01957
01958 else if( strcmp(punch.chPunch[ipPun],"LINP") == 0 )
01959 {
01960 if( strcmp(chTime,"LAST") != 0 )
01961 {
01962 static bool lgFirst=true;
01963
01964
01965
01966 PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[ipPun][0]);
01967 if( lgFirst )
01968 {
01969 lgFirst = false;
01970 PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[ipPun][0]);
01971 }
01972 }
01973 }
01974
01975 else if( strcmp(punch.chPunch[ipPun],"LINS") == 0 )
01976 {
01977
01978 if( strcmp(chTime,"LAST") != 0 )
01979 {
01980
01981 lgOK = true;
01982 punch_line(punch.ipPnunit[ipPun],"PUNS",lgOK,chDummy);
01983 }
01984 }
01985
01986 else if( strcmp(punch.chPunch[ipPun],"LINR") == 0 )
01987 {
01988
01989 if( strcmp(chTime,"LAST") != 0 )
01990 {
01991 Punch_Line_RT( punch.ipPnunit[ipPun] , "PUNC" );
01992 }
01993 }
01994
01995 else if( strcmp(punch.chPunch[ipPun],"LINA") == 0 )
01996 {
01997
01998 if( strcmp(chTime,"LAST") == 0 )
01999 {
02000
02001 for( j=0; j < LineSave.nsum; j++ )
02002 {
02003
02004
02005
02006 if( LineSv[j].xLineEnergy > 0. &&
02007 LineSv[j].sumlin[LineSave.lgLineEmergent] > 0. )
02008 {
02009
02010 fprintf( punch.ipPnunit[ipPun], "%12.5e",
02011 AnuUnit(LineSv[j].xLineEnergy) );
02012
02013 fprintf( punch.ipPnunit[ipPun], "\t%4.4s ",
02014 LineSv[j].chALab );
02015
02016 prt_wl( punch.ipPnunit[ipPun], LineSv[j].wavelength );
02017
02018 fprintf( punch.ipPnunit[ipPun], "\t%8.3f",
02019 log10(SDIV(LineSv[j].sumlin[0]) ) + radius.Conv2PrtInten );
02020
02021 fprintf( punch.ipPnunit[ipPun], "\t%8.3f",
02022 log10(SDIV(LineSv[j].sumlin[1]) ) + radius.Conv2PrtInten );
02023
02024 fprintf( punch.ipPnunit[ipPun], " \t%c\n",
02025 LineSv[j].chSumTyp);
02026 }
02027 }
02028 }
02029 }
02030
02031 else if( strcmp(punch.chPunch[ipPun],"LINI") == 0 )
02032 {
02033 if( strcmp(chTime,"LAST") == 0 &&
02034 (nzone/punch.LinEvery)*punch.LinEvery != nzone )
02035 {
02036
02037
02038 PunLineIntensity(punch.ipPnunit[ipPun]);
02039 }
02040 else if( strcmp(chTime,"LAST") != 0 )
02041 {
02042
02043 if( (punch.lgLinEvery && nzone == 1) ||
02044 (nzone/punch.LinEvery)*punch.LinEvery ==
02045 nzone )
02046 {
02047
02048
02049 PunLineIntensity(punch.ipPnunit[ipPun]);
02050 }
02051 }
02052 }
02053
02054 else if( strcmp( punch.chPunch[ipPun],"LEIL") == 0)
02055 {
02056
02057
02058 if( strcmp(chTime,"LAST") == 0 )
02059 {
02060 double absval , rel;
02061 long int n;
02062
02063
02064
02065
02066 const int NLINE_H2 = 31;
02067
02068 const int NLINE_NOTH_H2 = 5;
02069
02070 char chLabel[NLINE_NOTH_H2][5]=
02071 {"C 2", "O 1", "O 1","C 1", "C 1" };
02072 double Wl[NLINE_NOTH_H2]=
02073 {157.6 , 63.17 , 145.5 ,609.2 , 369.7 , };
02074
02075
02076
02077 double Wl_H2[NLINE_H2]=
02078 {2.121 ,
02079 28.21 , 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 ,
02080 4.408 , 4.180 , 3.996 , 3.845 , 3.724 , 3.626 , 3.547 , 3.485 , 3.437 , 3.403 ,
02081 3.381 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603};
02082
02083 for( n=0; n<NLINE_NOTH_H2; ++n )
02084 {
02085 fprintf(punch.ipPnunit[ipPun], "%s\t%.2f",chLabel[n] , Wl[n]);
02086
02087
02088 if( cdLine( chLabel[n] , (float)(Wl[n]*1e4) , &absval , &rel ) <= 0 )
02089 {
02090 fprintf(punch.ipPnunit[ipPun], " did not find\n");
02091 }
02092 else
02093 {
02094 fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
02095 }
02096 }
02097 fprintf(punch.ipPnunit[ipPun], "\n\n\n");
02098
02099
02100 if( h2.lgH2ON )
02101 {
02102 fprintf(punch.ipPnunit[ipPun],
02103 "Here are some of the H2 Intensities, The first one is the\n"
02104 "1-0 S(0) line and the following ones are the 0-0 S(X)\n"
02105 "lines where X goes from 0 to 29\n\n");
02106 for( n=0; n<NLINE_H2; ++n )
02107 {
02108 fprintf(punch.ipPnunit[ipPun], "%s\t%.3f","H2 " , Wl_H2[n]);
02109
02110 if( cdLine( "H2 " , (float)(Wl_H2[n]*1e4) , &absval , &rel ) <= 0 )
02111 {
02112 fprintf(punch.ipPnunit[ipPun], " did not find\n");
02113 }
02114 else
02115 {
02116 fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
02117 }
02118 }
02119 }
02120 }
02121 }
02122
02123 else if( strcmp( punch.chPunch[ipPun],"LEIS") == 0)
02124 {
02125 if( strcmp(chTime,"LAST") != 0 )
02126 {
02127
02128 double col_ci , col_oi , col_cii, col_heii;
02129 if( cdColm("carb" , 1 , &col_ci ) )
02130 TotalInsanity();
02131 if( cdColm("carb" , 2 , &col_cii ) )
02132 TotalInsanity();
02133 if( cdColm("oxyg" , 1 , &col_oi ) )
02134 TotalInsanity();
02135 if( cdColm("heli" , 2 , &col_heii ) )
02136 TotalInsanity();
02137
02138 fprintf( punch.ipPnunit[ipPun],
02139 "%.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"
02140 "%.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"
02141 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
02142 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
02143 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
02144
02145 radius.depth_mid_zone,
02146
02147 0.00,
02148
02149 rfield.extin_mag_V_point,
02150
02151 phycon.te ,
02152 dense.xIonDense[ipHYDROGEN][0],
02153 hmi.H2_total,
02154 dense.xIonDense[ipCARBON][0],
02155 dense.xIonDense[ipCARBON][1],
02156 dense.xIonDense[ipOXYGEN][0],
02157 findspecies("CO")->hevmol,
02158 findspecies("O2")->hevmol,
02159 findspecies("CH")->hevmol,
02160 findspecies("OH")->hevmol,
02161 dense.eden,
02162 dense.xIonDense[ipHELIUM][1],
02163 dense.xIonDense[ipHYDROGEN][1],
02164 hmi.Hmolec[ipMH3p],
02165 colden.colden[ipCOL_H0],
02166 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
02167 col_ci,
02168 col_cii,
02169 col_oi,
02170 findspecies("CO")->hevcol,
02171 findspecies("O2")->hevcol,
02172 findspecies("CH")->hevcol,
02173 findspecies("OH")->hevcol,
02174 colden.colden[ipCOL_elec],
02175 col_heii,
02176 colden.colden[ipCOL_Hp],
02177 colden.colden[ipCOL_H3p],
02178 hmi.H2_Solomon_dissoc_rate_used_H2g ,
02179 gv.rate_h2_form_grains_used_total,
02180 hmi.UV_Cont_rel2_Draine_DB96_depth,
02181
02182 CO_findrk("PHOTON,CO=>C,O"),
02183
02184 ionbal.PhotoRate_Shell[ipCARBON][0][2][0],
02185
02186 thermal.htot,
02187
02188 thermal.ctot,
02189
02190 thermal.heating[0][13],
02191
02192 MAX2(0.,gv.GasCoolColl),
02193
02194 -1.*MIN2(0.,gv.GasCoolColl),
02195
02196 thermal.heating[0][9],
02197
02198 hmi.HeatH2Dish_used,
02199
02200 hmi.HeatH2Dexc_used ,
02201
02202 thermal.heating[0][24] ,
02203
02204 thermal.heating[1][6] ,
02205
02206 thermal.heating[ipMAGNESIUM][0],
02207 thermal.heating[ipSULPHUR][0],
02208 thermal.heating[ipSILICON][0],
02209 thermal.heating[ipIRON][0],
02210 thermal.heating[ipSODIUM][0],
02211 thermal.heating[ipALUMINIUM][0],
02212 thermal.heating[ipCARBON][0],
02213 TauLines[ipT610].cool,
02214 TauLines[ipT370].cool,
02215 TauLines[ipT157].cool,
02216 TauLines[ipT63].cool,
02217 TauLines[ipT146].cool );
02218 }
02219 }
02220
02221 else if( strcmp( punch.chPunch[ipPun],"LLST") == 0)
02222 {
02223
02224
02225 if( strcmp(chTime,"LAST") == 0 )
02226 {
02227
02228 fprintf( punch.ipPnunit[ipPun],
02229 "%li" , iteration );
02230
02231
02232 if( punch.nLineList[ipPun] < 0 )
02233 TotalInsanity();
02234
02235 for( j=0; j<punch.nLineList[ipPun]; ++j )
02236 {
02237 double relative , absolute;
02238 if( (cdLine( punch.chLineListLabel[ipPun][j] ,
02239 punch.wlLineList[ipPun][j] ,
02240 &relative , &absolute ))<=0 )
02241 {
02242 fprintf(ioQQQ,"DISASTER - did not find a line in the Line List table\n");
02243 puts( "[Stop in PunchDo]" );
02244 cdEXIT(EXIT_FAILURE);
02245 }
02246
02247
02248
02249 if( punch.punarg[ipPun][0] > 0 )
02250 {
02251 fprintf( punch.ipPnunit[ipPun], "\t%.3f" , absolute );
02252 }
02253 else
02254 {
02255 fprintf( punch.ipPnunit[ipPun], "\t%.3e" , relative );
02256 }
02257 }
02258 fprintf( punch.ipPnunit[ipPun], "\n" );
02259 }
02260 }
02261
02262 else if( strcmp( punch.chPunch[ipPun],"RCO ") == 0)
02263 {
02264
02265 if( strcmp(chTime,"LAST") != 0 )
02266 {
02267 CO_punch_mol(punch.ipPnunit[ipPun],"CO",NULL,radius.depth_mid_zone);
02268 }
02269 }
02270
02271
02272 else if( strcmp( punch.chPunch[ipPun],"ROH ") == 0)
02273 {
02274
02275 if( strcmp(chTime,"LAST") != 0 )
02276 {
02277 CO_punch_mol(punch.ipPnunit[ipPun],"OH",NULL,radius.depth_mid_zone);
02278 }
02279 }
02280 else if( strcmp(punch.chPunch[ipPun],"MAP ") == 0 )
02281 {
02282
02283
02284 if( !map.lgMapDone &&
02285 (nzone == map.MapZone || strcmp(chTime,"LAST") == 0 ) )
02286 {
02287 lgTlkSav = called.lgTalk;
02288 called.lgTalk = true;
02289 map.lgMapBeingDone = true;
02290 map_do(punch.ipPnunit[ipPun] , " map");
02291 called.lgTalk = lgTlkSav;
02292 }
02293 }
02294
02295 else if( strcmp(punch.chPunch[ipPun],"MOLE") == 0 )
02296 {
02297 if( strcmp(chTime,"LAST") != 0 )
02298 {
02299
02300 fprintf( punch.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone );
02301
02302
02303 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
02304
02305
02306 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
02307
02308
02309 fprintf( punch.ipPnunit[ipPun], "%.5e\t" , CO_findrk("PHOTON,CO=>C,O") );
02310
02311
02312 fprintf( punch.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] );
02313
02314
02315 for(j=0; j<N_H_MOLEC; ++j )
02316 {
02317 fprintf(punch.ipPnunit[ipPun],"\t%.2e",hmi.Hmolec[j] );
02318 }
02319
02320
02321 for(j=0; j<mole.num_comole_calc; ++j )
02322 {
02323 fprintf(punch.ipPnunit[ipPun],"\t%.2e",COmole[j]->hevmol );
02324 }
02325 fprintf(punch.ipPnunit[ipPun],"\n");
02326 }
02327 }
02328
02329 else if( strcmp(punch.chPunch[ipPun],"OPAC") == 0 )
02330 {
02331
02332 if( strcmp(chTime,"LAST") == 0 )
02333 punch_opacity(punch.ipPnunit[ipPun],ipPun);
02334 }
02335
02336
02337 else if( strcmp(punch.chPunch[ipPun],"OPTc") == 0 )
02338 {
02339 if( strcmp(chTime,"LAST") == 0 )
02340 {
02341 for( j=0; j < rfield.nflux; j++ )
02342 {
02343 fprintf( punch.ipPnunit[ipPun],
02344 "%12.4e\t%.3e\t%12.4e\t%.3e\n",
02345 AnuUnit(rfield.AnuOrg[j]),
02346 opac.TauAbsFace[j]+opac.TauScatFace[j],
02347 opac.TauAbsFace[j],
02348 opac.TauScatFace[j] );
02349 }
02350 }
02351 }
02352
02353
02354 else if( strcmp(punch.chPunch[ipPun],"OPTf") == 0 )
02355 {
02356 if( strcmp(chTime,"LAST") == 0 )
02357 {
02358 long nu_hi , nskip;
02359 if( punch.punarg[ipPun][0] > 0. )
02360 {
02361
02362 j = ipFineCont( punch.punarg[ipPun][0] );
02363 }
02364 else
02365 {
02366 j = 0;
02367 }
02368
02369
02370 if( punch.punarg[ipPun][1]> 0. )
02371 {
02372 nu_hi = ipFineCont( punch.punarg[ipPun][1]);
02373 }
02374 else
02375 {
02376 nu_hi = rfield.nfine;
02377 }
02378 nskip = (long)punch.punarg[ipPun][2];
02379 do
02380 {
02381 float sum1=rfield.fine_opt_depth[j];
02382 float sum2=rfield.fine_opac_zone[j];
02383 float xnu = rfield.fine_anu[j+nskip/2];
02384 for( jj=1; jj<nskip; ++jj )
02385 {
02386 sum1 += rfield.fine_opt_depth[j+jj];
02387 sum2 += rfield.fine_opac_zone[j+jj];
02388 }
02389 if( sum2 > 0. )
02390 fprintf( punch.ipPnunit[ipPun],
02391 "%12.6e\t%.3e\t%.3e\n",
02392 AnuUnit(xnu),
02393 sum1/nskip ,
02394 sum2/nskip);
02395 j += nskip;
02396 }while( j < nu_hi );
02397 }
02398 }
02399
02400 else if( strcmp(punch.chPunch[ipPun]," OTS") == 0 )
02401 {
02402 ConMax = 0.;
02403 xLinMax = 0.;
02404 opConSum = 0.;
02405 opLinSum = 0.;
02406 ipLinMax = 1;
02407 ipConMax = 1;
02408
02409 for( j=0; j < rfield.nflux; j++ )
02410 {
02411 opConSum += rfield.otscon[j]*opac.opacity_abs[j];
02412 opLinSum += rfield.otslin[j]*opac.opacity_abs[j];
02413 if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax )
02414 {
02415 xLinMax = rfield.otslin[j]*opac.opacity_abs[j];
02416 ipLinMax = j+1;
02417 }
02418 if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax )
02419 {
02420 ConMax = rfield.otscon[j]*opac.opacity_abs[j];
02421 ipConMax = j+1;
02422 }
02423 }
02424 fprintf( punch.ipPnunit[ipPun],
02425 "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n",
02426 opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1]
02427 , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1]
02428 , rfield.anu[ipConMax-1], ConMax );
02429 }
02430
02431 else if( strcmp(punch.chPunch[ipPun],"OVER") == 0 )
02432 {
02433
02434
02435 double toosmall = SMALLFLOAT ,
02436 hold;
02437
02438
02439
02440 if( strcmp(chTime,"LAST") != 0 )
02441 {
02442
02443
02444 fprintf( punch.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone );
02445
02446
02447 if(dynamics.Cool > dynamics.Heat)
02448 {
02449 fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f",
02450 log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat)) );
02451 }
02452 else
02453 {
02454 double diff = fabs(thermal.htot-dynamics.Cool);
02455 fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f",
02456 log10(phycon.te), log10(SDIV(diff)) );
02457 }
02458
02459
02460 fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f",
02461 log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) );
02462
02463
02464 fprintf( punch.ipPnunit[ipPun], "\t%.4f",
02465
02466 log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) );
02467
02468
02469 fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f",
02470 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])),
02471 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );
02472
02473
02474 for( j=1; j <= 3; j++ )
02475 {
02476 double arg1 = SDIV(dense.gas_phase[ipHELIUM]);
02477 arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 );
02478 fprintf( punch.ipPnunit[ipPun], "\t%.4f",
02479 log10(arg1) );
02480 }
02481
02482
02483 hold = SDIV(dense.gas_phase[ipCARBON]);
02484 hold = findspecies("CO")->hevmol/hold;
02485 hold = MAX2(toosmall, hold );
02486 fprintf( punch.ipPnunit[ipPun], "\t%.4f", log10(hold) );
02487
02488
02489 for( j=1; j <= 4; j++ )
02490 {
02491 hold = SDIV(dense.gas_phase[ipCARBON]);
02492 hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold);
02493 fprintf( punch.ipPnunit[ipPun], "\t%.4f",
02494 log10(hold) );
02495 }
02496
02497
02498 for( j=1; j <= 6; j++ )
02499 {
02500 hold = SDIV(dense.gas_phase[ipOXYGEN]);
02501 hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold);
02502 fprintf( punch.ipPnunit[ipPun], "\t%.4f",
02503 log10(hold) );
02504 }
02505
02506
02507 fprintf( punch.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point);
02508
02509
02510 fprintf( punch.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended);
02511 }
02512 }
02513
02514 else if( strcmp(punch.chPunch[ipPun]," PDR") == 0 )
02515 {
02516
02517 if( strcmp(chTime,"LAST") != 0 )
02518 {
02519
02520
02521
02522
02523
02524
02525
02526 fprintf( punch.ipPnunit[ipPun],
02527 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t",
02528 radius.depth_mid_zone,
02529
02530 colden.colden[ipCOL_HTOT],
02531 phycon.te,
02532
02533 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
02534
02535 2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN],
02536 2.*hmi.Hmolec[ipMH2s]/dense.gas_phase[ipHYDROGEN],
02537
02538 dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]),
02539 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
02540 findspecies("H2O")->hevmol/SDIV(dense.gas_phase[ipOXYGEN]),
02541
02542 hmi.UV_Cont_rel2_Habing_TH85_depth);
02543
02544
02545 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
02546
02547
02548 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
02549
02550
02551 fprintf( punch.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] );
02552 }
02553 }
02554
02555 else if( strcmp(punch.chPunch[ipPun],"PHYS") == 0 )
02556 {
02557 if( strcmp(chTime,"LAST") != 0 )
02558 {
02559
02560 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
02561 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN],
02562 dense.eden, thermal.htot, wind.AccelTot, geometry.FillFac );
02563 }
02564 }
02565
02566 else if( strcmp(punch.chPunch[ipPun],"PRES") == 0 )
02567 {
02568
02569 if( strcmp(chTime,"LAST") != 0 )
02570 {
02571 fprintf( punch.ipPnunit[ipPun],
02572 "%.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",
02573
02574 radius.depth_mid_zone,
02575
02576 pressure.PresTotlCorrect,
02577
02578 pressure.PresTotlCurr,
02579
02580
02581
02582 pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon,
02583
02584 pressure.PresTotlInit,
02585
02586 pressure.PresGasCurr,
02587
02588 pressure.PresRamCurr,
02589
02590 pressure.pres_radiation_lines_curr,
02591
02592 pressure.PresInteg - pressure.pinzon,
02593
02594 wind.windv/1e5,
02595
02596 timesc.sound_speed_adiabatic/1e5,
02597
02598 magnetic.pressure ,
02599
02600 DoppVel.TurbVel/1e5 ,
02601
02602 pressure.PresTurbCurr*DoppVel.lgTurb_pressure ,
02603 TorF(conv.lgConvPres) );
02604 }
02605 }
02606
02607 else if( punch.chPunch[ipPun][0]=='R' )
02608 {
02609
02610 if( strcmp(punch.chPunch[ipPun],"RADI") == 0 )
02611 {
02612
02613 if( strcmp(chTime,"LAST") != 0 )
02614 {
02615 fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
02616 nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
02617 radius.drad );
02618 }
02619 }
02620
02621 else if( strcmp(punch.chPunch[ipPun],"RADO") == 0 )
02622 {
02623
02624 if( strcmp(chTime,"LAST") == 0 )
02625 {
02626 fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
02627 nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
02628 radius.drad );
02629 }
02630 }
02631
02632 else if( strcmp(punch.chPunch[ipPun],"RESU") == 0 )
02633 {
02634
02635 if( strcmp(chTime,"LAST") == 0 )
02636 punResults(punch.ipPnunit[ipPun]);
02637 }
02638 else
02639 {
02640
02641 TotalInsanity();
02642 }
02643 }
02644
02645 else if( strcmp(punch.chPunch[ipPun],"SECO") == 0 )
02646 {
02647
02648 if( strcmp(chTime,"LAST") != 0 )
02649 fprintf(punch.ipPnunit[ipPun],
02650 "%.5e\t%.3e\t%.3e\t%.3e\n",
02651 radius.depth ,
02652 secondaries.csupra[ipHYDROGEN][0],
02653 secondaries.csupra[ipHYDROGEN][0]*2.02,
02654 secondaries.x12tot );
02655 }
02656
02657 else if( strcmp(punch.chPunch[ipPun],"SOUS") == 0 )
02658 {
02659
02660
02661 if( strcmp(chTime,"LAST") != 0 )
02662 {
02663 limit = MIN2(rfield.ipMaxBolt,rfield.nflux);
02664 for( j=0; j < limit; j++ )
02665 {
02666 fprintf( punch.ipPnunit[ipPun],
02667 "%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
02668 rfield.anu[j],
02669 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j],
02670 opac.opacity_abs[j],
02671 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j]),
02672 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j])/plankf(j) );
02673 }
02674 }
02675 }
02676
02677 else if( strcmp(punch.chPunch[ipPun],"SOUD") == 0 )
02678 {
02679
02680
02681 j = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] + 2;
02682 fprintf( punch.ipPnunit[ipPun],
02683 "%.4e\t%.4e\t%.4e\t%.4e\n",
02684 opac.TauAbsFace[j-1],
02685 rfield.ConEmitLocal[nzone][j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]),
02686 rfield.otscon[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1],
02687 rfield.otscon[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/opac.opacity_abs[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] );
02688 }
02689
02690
02691 else if( strcmp(punch.chPunch[ipPun],"SPEC") == 0 )
02692 {
02693 PunchSpecial(punch.ipPnunit[ipPun],chTime);
02694 }
02695
02696 else if( strcmp(punch.chPunch[ipPun],"TEGR") == 0 )
02697 {
02698
02699 if( strcmp(chTime,"LAST") != 0 )
02700 {
02701 fprintf( punch.ipPnunit[ipPun], " " );
02702 for( j=0; j < NGRID; j++ )
02703 {
02704 fprintf( punch.ipPnunit[ipPun],
02705 "%4ld%11.3e%11.3e%11.3e\n",
02706 thermal.nZonGrid[j],
02707 thermal.TeGrid[j],
02708 thermal.HtGrid[j],
02709 thermal.ClGrid[j] );
02710 }
02711 }
02712 }
02713
02714 else if( strcmp(punch.chPunch[ipPun],"TEMP") == 0 )
02715 {
02716 static double deriv_old=-1;
02717 double deriv=-1. , deriv_sec;
02718
02719
02720 if( nzone >1 )
02721 {
02722 deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
02723 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e",
02724 radius.depth_mid_zone,
02725 phycon.te,
02726 deriv );
02727
02728 if( nzone > 2 )
02729 {
02730 deriv_sec = (deriv-deriv_old)/ radius.drad;
02731 fprintf( punch.ipPnunit[ipPun], "\t%.2e",
02732 deriv_sec );
02733 }
02734 deriv_old = deriv;
02735 fprintf( punch.ipPnunit[ipPun], "\n");
02736 }
02737 }
02738
02739
02740 else if( strcmp(punch.chPunch[ipPun],"TIMD") == 0 )
02741 {
02742 if( strcmp(chTime,"LAST") == 0 )
02743 DynaPunchTimeDep( punch.ipPnunit[ipPun] , "END" );
02744 }
02745
02746 else if( strcmp(punch.chPunch[ipPun],"XTIM") == 0 )
02747 {
02748 static double ElapsedTime , ZoneTime;
02749 if( nzone<=1 )
02750 {
02751 ElapsedTime = cdExecTime();
02752 ZoneTime = 0.;
02753 }
02754 else
02755 {
02756 double t = cdExecTime();
02757 ZoneTime = t - ElapsedTime;
02758 ElapsedTime = t;
02759 }
02760
02761
02762 fprintf( punch.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n",
02763 nzone, ZoneTime , ElapsedTime );
02764 }
02765
02766 else if( strcmp(punch.chPunch[ipPun],"TPRE") == 0 )
02767 {
02768
02769 fprintf( punch.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n",
02770 nzone, phycon.TeInit, phycon.TeProp, phycon.te,
02771 (phycon.TeProp- phycon.te)/phycon.te );
02772 }
02773
02774 else if( strcmp(punch.chPunch[ipPun],"WIND") == 0 )
02775 {
02776
02777
02778 if( punch.punarg[ipPun][0] == 0 &&
02779 strcmp(chTime,"LAST") == 0 )
02780 {
02781
02782 fprintf( punch.ipPnunit[ipPun],
02783 "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
02784 radius.Radius_mid_zone,
02785 radius.depth_mid_zone,
02786 wind.windv,
02787 wind.AccelTot,
02788 wind.AccelLine,
02789 wind.AccelCont ,
02790 wind.fmul );
02791 }
02792 else if( punch.punarg[ipPun][0] == 1 &&
02793 strcmp(chTime,"LAST") != 0 )
02794 {
02795
02796 fprintf( punch.ipPnunit[ipPun],
02797 "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
02798 radius.Radius_mid_zone,
02799 radius.depth_mid_zone,
02800 wind.windv,
02801 wind.AccelTot,
02802 wind.AccelLine,
02803 wind.AccelCont ,
02804 wind.fmul );
02805 }
02806 }
02807 else if( strcmp(punch.chPunch[ipPun],"XATT") == 0 )
02808 {
02809
02810 ASSERT( grid.lgOutputTypeOn[2] );
02811
02812 if( strcmp(chTime,"LAST") == 0 )
02813 {
02814 if( grid.lgGrid )
02815 punchFITSfile( punch.ipPnunit[ipPun], 2 );
02816 else
02817 {
02818 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02819 puts( "[Stop in PunchDo]" );
02820 cdEXIT(EXIT_FAILURE);
02821 }
02822 }
02823 }
02824 else if( strcmp(punch.chPunch[ipPun],"XRFI") == 0 )
02825 {
02826
02827 ASSERT( grid.lgOutputTypeOn[3] );
02828
02829 if( strcmp(chTime,"LAST") == 0 )
02830 {
02831 if( grid.lgGrid )
02832 punchFITSfile( punch.ipPnunit[ipPun], 3 );
02833 else
02834 {
02835 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02836 puts( "[Stop in PunchDo]" );
02837 cdEXIT(EXIT_FAILURE);
02838 }
02839 }
02840 }
02841 else if( strcmp(punch.chPunch[ipPun],"XINC") == 0 )
02842 {
02843
02844 ASSERT( grid.lgOutputTypeOn[1] );
02845
02846 if( strcmp(chTime,"LAST") == 0 )
02847 {
02848 if( grid.lgGrid )
02849 punchFITSfile( punch.ipPnunit[ipPun], 1 );
02850 else
02851 {
02852 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02853 puts( "[Stop in PunchDo]" );
02854 cdEXIT(EXIT_FAILURE);
02855 }
02856 }
02857 }
02858 else if( strcmp(punch.chPunch[ipPun],"XDFR") == 0 )
02859 {
02860
02861 ASSERT( grid.lgOutputTypeOn[5] );
02862
02863 if( strcmp(chTime,"LAST") == 0 )
02864 {
02865 if( grid.lgGrid )
02866 punchFITSfile( punch.ipPnunit[ipPun], 5 );
02867 else
02868 {
02869 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02870 puts( "[Stop in PunchDo]" );
02871 cdEXIT(EXIT_FAILURE);
02872 }
02873 }
02874 }
02875 else if( strcmp(punch.chPunch[ipPun],"XDFO") == 0 )
02876 {
02877
02878 ASSERT( grid.lgOutputTypeOn[4] );
02879
02880 if( strcmp(chTime,"LAST") == 0 )
02881 {
02882 if( grid.lgGrid )
02883 punchFITSfile( punch.ipPnunit[ipPun], 4 );
02884 else
02885 {
02886 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02887 puts( "[Stop in PunchDo]" );
02888 cdEXIT(EXIT_FAILURE);
02889 }
02890 }
02891 }
02892 else if( strcmp(punch.chPunch[ipPun],"XLNR") == 0 )
02893 {
02894
02895 ASSERT( grid.lgOutputTypeOn[7] );
02896
02897 if( strcmp(chTime,"LAST") == 0 )
02898 {
02899 if( grid.lgGrid )
02900 punchFITSfile( punch.ipPnunit[ipPun], 7 );
02901 else
02902 {
02903 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02904 puts( "[Stop in PunchDo]" );
02905 cdEXIT(EXIT_FAILURE);
02906 }
02907 }
02908 }
02909 else if( strcmp(punch.chPunch[ipPun],"XLNO") == 0 )
02910 {
02911
02912 ASSERT( grid.lgOutputTypeOn[6] );
02913
02914 if( strcmp(chTime,"LAST") == 0 )
02915 {
02916 if( grid.lgGrid )
02917 punchFITSfile( punch.ipPnunit[ipPun], 6 );
02918 else
02919 {
02920 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02921 puts( "[Stop in PunchDo]" );
02922 cdEXIT(EXIT_FAILURE);
02923 }
02924 }
02925 }
02926 else if( strcmp(punch.chPunch[ipPun],"XREF") == 0 )
02927 {
02928
02929 ASSERT( grid.lgOutputTypeOn[9] );
02930
02931 if( strcmp(chTime,"LAST") == 0 )
02932 {
02933 if( grid.lgGrid )
02934 punchFITSfile( punch.ipPnunit[ipPun], 9 );
02935 else
02936 {
02937 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02938 puts( "[Stop in PunchDo]" );
02939 cdEXIT(EXIT_FAILURE);
02940 }
02941 }
02942 }
02943 else if( strcmp(punch.chPunch[ipPun],"XTRN") == 0 )
02944 {
02945
02946 ASSERT( grid.lgOutputTypeOn[8] );
02947
02948 if( strcmp(chTime,"LAST") == 0 )
02949 {
02950 if( grid.lgGrid )
02951 punchFITSfile( punch.ipPnunit[ipPun], 8 );
02952 else
02953 {
02954 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02955 puts( "[Stop in PunchDo]" );
02956 cdEXIT(EXIT_FAILURE);
02957 }
02958 }
02959 }
02960 else if( strcmp(punch.chPunch[ipPun],"XSPM") == 0 )
02961 {
02962
02963 ASSERT( grid.lgOutputTypeOn[10] );
02964
02965 if( strcmp(chTime,"LAST") == 0 )
02966 {
02967 if( grid.lgGrid )
02968 punchFITSfile( punch.ipPnunit[ipPun], 10 );
02969 else
02970 {
02971 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
02972 puts( "[Stop in PunchDo]" );
02973 cdEXIT(EXIT_FAILURE);
02974 }
02975 }
02976 }
02977 else
02978 {
02979
02980 fprintf( ioQQQ, " PunchDo does not recognize flag %4.4s set by ParsePunch. This is impossible.\n",
02981 punch.chPunch[ipPun] );
02982 ShowMe();
02983 puts( "[Stop in PunchDo]" );
02984 cdEXIT(EXIT_FAILURE);
02985 }
02986
02987
02988
02989
02990 if( strcmp(chTime,"LAST") == 0 && !iterations.lgLastIt && punch.lgHashEndIter &&
02991
02992
02993 punch.lg_separate_iterations[ipPun] )
02994 {
02995 if( dynamics.lgStatic && strcmp( punch.chHashString , "TIME_DEP" )==0 )
02996 {
02997 fprintf( punch.ipPnunit[ipPun], "\"time=%f\n",
02998 dynamics.time_elapsed );
02999 }
03000 else
03001 {
03002 fprintf( punch.ipPnunit[ipPun], "%s\n",
03003 punch.chHashString );
03004 }
03005 }
03006 if( punch.lgFLUSH )
03007 fflush( punch.ipPnunit[ipPun] );
03008 }
03009 }
03010
03011 DEBUG_EXIT( "PunchDo()" );
03012 return;
03013 }
03014
03015
03016 static void PunLineIntensity(FILE * ioPUN)
03017 {
03018 char chCard[INPUT_LINE_LENGTH];
03019 bool lgEOF;
03020 long int i;
03021
03022 DEBUG_ENTRY( "PunLineIntensity()" );
03023
03024
03025
03026
03027
03028 input_init();
03029 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
03030
03031 lgEOF = false;
03032 while( !lgEOF )
03033 {
03034 input_readarray(chCard,&lgEOF);
03035 if( !lgEOF )
03036 {
03037 fprintf( ioPUN, "%s\n", chCard );
03038 }
03039 }
03040
03041
03042 cdWarnings( ioPUN);
03043 cdCautions( ioPUN);
03044 fprintf( ioPUN, "zone=%5ld\n", nzone );
03045 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
03046 fprintf( ioPUN, "begin emission lines\n" );
03047
03048
03049 PunResults1Line(ioPUN," ",0,0.,"Start");
03050 for( i=0; i < LineSave.nsum; i++ )
03051 {
03052 if( LineSv[i].sumlin[LineSave.lgLineEmergent] > 0. )
03053 {
03054 PunResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
03055 LineSv[i].sumlin[LineSave.lgLineEmergent], "Line ");
03056 }
03057 }
03058
03059 PunResults1Line(ioPUN," ",0,0.,"Flush");
03060
03061 fprintf( ioPUN, " \n" );
03062 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
03063
03064
03065 DEBUG_EXIT( "PunLineIntensity()" );
03066 return;
03067 }
03068
03069
03070 static bool lgPopsFirstCall , lgPunchOpticalDepths;
03071
03072
03073 static void PunchLineStuff(
03074 FILE * ioPUN,
03075 const char *chJob ,
03076 float xLimit )
03077 {
03078
03079 long int nelem , ipLo , ipHi , i , ipISO;
03080 long index = 0;
03081 static bool lgFirst=true;
03082
03083 DEBUG_ENTRY( "PunchLineStuff()" );
03084
03085
03086 if( strcmp( &*chJob , "optical" ) == 0 )
03087 {
03088
03089 lgPunchOpticalDepths = true;
03090 lgPopsFirstCall = false;
03091 }
03092 else if( strcmp( &*chJob , "populat" ) == 0 )
03093 {
03094 lgPunchOpticalDepths = false;
03095
03096 if( lgFirst )
03097 {
03098 lgPopsFirstCall = true;
03099 fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
03100 lgFirst = false;
03101 }
03102 else
03103 {
03104 lgPopsFirstCall = false;
03105 }
03106 }
03107 else
03108 {
03109 fprintf( ioQQQ, " insane job in PunchLineStuff =%s\n",
03110 &*chJob );
03111 puts( "[Stop in PunchLineStuff]" );
03112 cdEXIT(EXIT_FAILURE);
03113 }
03114
03115
03116
03117
03118 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
03119 {
03120 for( nelem=ipISO; nelem < LIMELM; nelem++ )
03121 {
03122 if( dense.lgElmtOn[nelem] )
03123 {
03124
03125 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
03126 {
03127 for( ipLo=0; ipLo <ipHi; ipLo++ )
03128 {
03129 ++index;
03130 pun1Line( &EmisLines[ipISO][nelem][ipHi][ipLo] , ioPUN , xLimit , index , dense.xIonDense[nelem][nelem+1-ipISO] );
03131 }
03132 }
03133
03134
03135
03136 if( lgPunchOpticalDepths )
03137 {
03138
03139
03140
03141
03142
03143 for( ipHi=iso.quant_desig[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1].n+1; ipHi < iso.nLyman[ipISO]; ipHi++ )
03144 {
03145 ++index;
03146 pun1Line( &iso.ExtraLymanLines[ipISO][nelem][ipHi] , ioPUN , xLimit , index , dense.xIonDense[nelem][nelem+1-ipISO]);
03147 }
03148 }
03149 }
03150 }
03151 }
03152
03153
03154 for( i=1; i < nLevel1; i++ )
03155 {
03156 ++index;
03157 pun1Line( &TauLines[i] , ioPUN , xLimit , index , 1. );
03158 }
03159
03160 for( i=0; i < nWindLine; i++ )
03161 {
03162 if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
03163 {
03164 ++index;
03165 pun1Line( &TauLine2[i] , ioPUN , xLimit , index , 1. );
03166 }
03167 }
03168
03169 for( i=0; i < nUTA; i++ )
03170 {
03171 ++index;
03172 pun1Line( &UTALines[i] , ioPUN , xLimit , index , 1. );
03173 }
03174
03175
03176
03177 FeIIPunchLineStuff( ioPUN , xLimit , index);
03178
03179
03180 H2_PunchLineStuff( ioPUN , xLimit , index);
03181
03182
03183 for( i=0; i < nCORotate; i++ )
03184 {
03185 ++index;
03186 pun1Line( &C12O16Rotate[i] , ioPUN , xLimit , index , 1. );
03187 }
03188
03189
03190 for( i=0; i < nCORotate; i++ )
03191 {
03192 ++index;
03193 pun1Line( &C13O16Rotate[i] , ioPUN , xLimit , index , 1. );
03194 }
03195
03196 fprintf( ioPUN , "%s\n",punch.chHashString );
03197
03198 DEBUG_EXIT( "PunchLineStuff()" );
03199 return;
03200 }
03201
03202
03203 void pun1Line( EmLine * t , FILE * ioPUN , float xLimit , long index , double abundance)
03204 {
03205
03206 if( lgPunchOpticalDepths )
03207 {
03208
03209 if( t->TauIn >= xLimit )
03210 {
03211
03212 fprintf( ioPUN , "%2.2s%2.2s\t",
03213 elementnames.chElementSym[t->nelem-1] ,
03214 elementnames.chIonStage[t->IonStg-1] );
03215
03216
03217
03218 if( strcmp( punch.chConPunEnr[punch.ipConPun], "labl" )== 0 )
03219 {
03220 prt_wl( ioPUN , t->WLAng );
03221 }
03222 else
03223 {
03224
03225 fprintf( ioPUN , "%.7e", AnuUnit((float)(t->EnergyWN * WAVNRYD)) );
03226 }
03227
03228 fprintf( ioPUN , "\t%.3f", t->TauIn );
03229
03230 fprintf(ioPUN, "\t%.3e",
03231 t->dampXvel / DoppVel.doppler[ t->nelem -1 ] );
03232 fprintf(ioPUN, "\n");
03233 }
03234 }
03235 else if( lgPopsFirstCall )
03236 {
03237 char chLbl[11];
03238
03239 strcpy( chLbl, chLineLbl(t) );
03240 fprintf(ioPUN, "%li\t%s" , index , chLbl );
03241
03242 fprintf(ioPUN, "\t%.0f\t%.0f",
03243 t->gLo ,t->gHi);
03244
03245 fprintf(ioPUN, "\t%.2f\t%.3e",
03246 t->EnergyWN ,t->gf);
03247 fprintf(ioPUN, "\n");
03248 }
03249 else
03250 {
03251
03252 if( t->PopHi > xLimit )
03253 {
03254 fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index,
03255
03256
03257
03258
03259 t->PopLo*abundance ,
03260 t->PopHi*abundance );
03261 }
03262 }
03263 }
03264
03265
03266 static void PunchNewContinuum(FILE * ioPUN, long ipCon )
03267 {
03268 long int ipLo, ipHi,
03269 j ,
03270 ncells;
03271
03272 double
03273 wllo=3500. ,
03274 wlhi=7000. ,
03275 resolution = 2.;
03276
03277 double *energy,
03278 *cont_incid,
03279 *cont_atten,
03280 *diffuse_in,
03281 *diffuse_out,
03282 *emis_lines_out,
03283 *emis_lines_in;
03284
03285
03286 wllo = MAX2( rfield.anu[0] , punch.cp_range_min[ipCon] );
03287 if( punch.cp_range_max[ipCon] > 0. )
03288 {
03289
03290 wlhi = MIN2( rfield.anu[rfield.nflux-1] , punch.cp_range_max[ipCon] );
03291 }
03292 else
03293 {
03294
03295 wlhi = rfield.anu[rfield.nflux-1];
03296 }
03297
03298 if( punch.cp_resolving_power[ipCon] != 0. )
03299 {
03300
03301 ipLo = LONG_MAX;
03302 ipHi = LONG_MAX;
03303
03304 ncells = (long)((wlhi-wllo)/resolution + 1);
03305 }
03306 else
03307 {
03308
03309 ipLo = ipoint(wllo)-1;
03310 ipHi = ipoint(wlhi)-1;
03311 ncells = ipHi - ipLo + 1;
03312 }
03313
03314
03315 energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
03316 cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
03317 cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
03318 diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
03319 diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
03320 emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
03321 emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
03322
03323
03324
03325
03326 if( punch.cp_resolving_power[ipCon] != 0. )
03327 {
03328
03329 energy[0] = wlhi;
03330 j = 0;
03331 while( energy[j]-resolution > wllo )
03332 {
03333 ++j;
03334 ASSERT( j< ncells+1 );
03335 energy[j] = energy[j-1] - resolution;
03336 }
03337
03338 ncells = j;
03339
03340 for( j=0; j<ncells; ++j )
03341 {
03342 energy[j] = RYDLAM / energy[j];
03343 }
03344 }
03345 else
03346 {
03347
03348 for( j=0; j<ncells; ++j )
03349 {
03350 energy[j] = rfield.AnuOrg[j+ipLo] - rfield.widflx[j+ipLo]/2.;
03351 }
03352 }
03353
03354
03355
03356 cdSPEC( 1 , energy , ncells , cont_incid );
03357
03358 cdSPEC( 2 , energy , ncells , cont_atten );
03359
03360 cdSPEC( 5 , energy , ncells , diffuse_in );
03361
03362 cdSPEC( 4 , energy , ncells , diffuse_out );
03364
03365 cdSPEC( 6 , energy , ncells , emis_lines_out );
03366 cdSPEC( 7 , energy , ncells , emis_lines_in );
03367
03368
03369
03370
03371 for( j=0; j<ncells-1; ++j )
03372 {
03373
03374 fprintf( ioPUN,"%.3e\t", AnuUnit((float)(energy[j]+rfield.widflx[j+ipLo]/2.) ) );
03375 fprintf( ioPUN,"%.3e\t", cont_incid[j] );
03376 fprintf( ioPUN,"%.3e\t", cont_atten[j] );
03377 fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
03378 fprintf( ioPUN,"%.3e",
03379 emis_lines_out[j]+emis_lines_in[j] );
03380 fprintf( ioPUN,"\n" );
03381 }
03382
03383 free(energy);
03384 free(cont_incid);
03385 free(diffuse_in);
03386 free(diffuse_out);
03387 free(cont_atten);
03388 free(emis_lines_out);
03389 free(emis_lines_in);
03390
03391
03392 }
03393
03394
03395 static void AGN_Hemis(FILE *ioPUN )
03396 {
03397 const int NTE = 4;
03398 float te[NTE] = {5000., 10000., 15000., 20000.};
03399 float *agn_continuum[NTE];
03400 double TempSave = phycon.te;
03401 long i , j;
03402
03403
03404
03405 for( i=0;i<NTE; ++i)
03406 {
03407 agn_continuum[i] = (float *)MALLOC((unsigned)rfield.nflux*sizeof(float) );
03408
03409
03410 phycon.te = te[i];
03411
03412 tfidle(true);
03413
03414
03415 ConvPresTempEdenIoniz();
03416
03417
03418
03419 RT_diffuse();
03420 for(j=0;j<rfield.nflux; ++j )
03421 {
03422 agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(float)dense.eden/
03423 (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] );
03424 }
03425 }
03426
03427
03428 fprintf(ioPUN,"wl");
03429 for( i=0;i<NTE; ++i)
03430 {
03431 fprintf(ioPUN,"\tT=%.0f",te[i]);
03432 }
03433 fprintf( ioPUN , "\tcont\n");
03434
03435
03436 for(j=0;j<rfield.nflux; ++j )
03437 {
03438 fprintf( ioPUN , "%12.5e",
03439 AnuUnit(rfield.AnuOrg[j]) );
03440
03441 for( i=0;i<NTE; ++i)
03442 {
03443 fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
03444 }
03445
03446 fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j]);
03447 }
03448
03449
03450 for( i=0;i<NTE; ++i)
03451 {
03452 free( agn_continuum[i] );
03453 }
03454
03455
03456 phycon.te = TempSave;
03457
03458 tfidle(true);
03459
03460 fprintf( ioQQQ, "AGN_Hemis - result of punch AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
03461 puts( "[Stop in PunchDo]" );
03462 cdEXIT(EXIT_FAILURE);
03463 }
03464
03465
03466
03467 static void punResults(FILE* ioPUN)
03468 {
03469 char chCard[INPUT_LINE_LENGTH];
03470 bool lgEOF;
03471 long int i , nelem , ion;
03472
03473 DEBUG_ENTRY( "punResults()" );
03474
03475
03476
03477 lgEOF = false;
03478
03479 input_init();
03480
03481 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
03482 lgEOF = false;
03483 input_readarray(chCard,&lgEOF);
03484 while( !lgEOF )
03485 {
03486
03487 char chCAPS[INPUT_LINE_LENGTH];
03488 strcpy( chCAPS , chCard );
03489 caps( chCAPS );
03490
03491 if( !nMatch( "HIDE" , chCAPS ) )
03492 fprintf( ioPUN, "%s\n", chCard );
03493 input_readarray(chCard,&lgEOF);
03494 }
03495
03496
03497 cdWarnings(ioPUN);
03498 cdCautions(ioPUN);
03499 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
03500
03501 fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );
03502
03503 fprintf( ioPUN, "BEGIN EMISSION LINES\n" );
03504 PunResults1Line(ioPUN," ",0,0.,"Start");
03505
03506 for( i=0; i < LineSave.nsum; i++ )
03507 {
03508 if( LineSv[i].sumlin[LineSave.lgLineEmergent] > 0. )
03509 {
03510 PunResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
03511 LineSv[i].sumlin[LineSave.lgLineEmergent],"Line ");
03512 }
03513 }
03514
03515 PunResults1Line(ioPUN," ",0,0.,"Flush");
03516
03517 fprintf( ioPUN, " \n" );
03518
03519 fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" );
03520
03521
03522
03523
03524 ASSERT( LIMELM == 30 );
03525
03526
03527 for( nelem=0; nelem<LIMELM; nelem++ )
03528 {
03529 for(ion=0; ion < nelem+1; ion++)
03530 {
03531 fprintf( ioPUN, " %10.3e", mean.xIonMeans[0][nelem][ion] );
03532
03533 if( nelem==9|| nelem==19 || nelem==29 )
03534 {
03535 fprintf( ioPUN, "\n" );
03536 }
03537 }
03538 }
03539
03540 fprintf( ioPUN, "END OF RESULTS\n" );
03541 fprintf( ioPUN, "**********************************************************************************************************************************\n" );
03542
03543 DEBUG_EXIT( "punResults()" );
03544 return;
03545 }
03546
03547 static const int LINEWIDTH = 6;
03548
03549
03550
03551 static void PunResults1Line(
03552 FILE * ioPUN,
03553
03554 const char *chLab,
03555 float wl,
03556 double xInten,
03557
03558 const char *chFunction)
03559 {
03560
03561 long int i;
03562 static float wavelength[LINEWIDTH];
03563 static long ipLine;
03564 static double xIntenSave[LINEWIDTH];
03565 static char chLabSave[LINEWIDTH][5];
03566
03567 DEBUG_ENTRY( "PunResults1Line()" );
03568
03569
03570
03571 if( strcmp(chFunction,"Start") == 0 )
03572 {
03573 ipLine = 0;
03574 }
03575 else if( strcmp(chFunction,"Line ") == 0 )
03576 {
03577
03578 wavelength[ipLine] = wl;
03579 strcpy( chLabSave[ipLine], chLab );
03580 xIntenSave[ipLine] = xInten;
03581
03582
03583
03584 ++ipLine;
03585
03586
03587 if( ( strcmp(punch.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
03588 {
03589
03590 for( i=0; i < ipLine; i++ )
03591 {
03592 fprintf( ioPUN, " %4.4s ", chLabSave[i] );
03593 prt_wl( ioPUN, wavelength[i] );
03594 fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
03595
03596
03597 if( strcmp(punch.chPunRltType,"column") == 0 )
03598 fprintf( ioPUN, "\n" );
03599 }
03600
03601 if( strcmp(punch.chPunRltType,"array ") == 0 )
03602 fprintf( ioPUN, " \n" );
03603 ipLine = 0;
03604 }
03605 }
03606 else if( strcmp(chFunction,"Flush") == 0 )
03607 {
03608 if( ipLine > 0 )
03609 {
03610
03611
03612
03613
03614
03615 for( i=0; i < ipLine; i++ )
03616 {
03617 fprintf( ioPUN, " %4.4s", chLabSave[i] );
03618 prt_wl( ioPUN, wavelength[i] );
03619 fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
03620
03621
03622 if( strcmp(punch.chPunRltType,"column") == 0 )
03623 fprintf( ioPUN, "\n" );
03624 }
03625 if( strcmp(punch.chPunRltType,"array ") == 0 )
03626 fprintf( ioPUN, " \n" );
03627 }
03628 }
03629 else
03630 {
03631 fprintf( ioQQQ, " PunResults1Line called with insane request=%5.5s\n",
03632 chFunction );
03633 puts( "[Stop in result]" );
03634 cdEXIT(EXIT_FAILURE);
03635 }
03636
03637 DEBUG_EXIT( "PunResults1Line()" );
03638 return;
03639 }
03640
03641 static const int NENR_GAUNT = 37;
03642 static const int NTE_GAUNT = 21;
03643
03644
03645 static void PunchGaunts(FILE* ioPUN)
03646 {
03647 long int i,
03648 charge,
03649 ite,
03650 j;
03651
03652 float ener[NENR_GAUNT],
03653 ste[NTE_GAUNT],
03654 z;
03655 double tempsave;
03656 double g[NENR_GAUNT][NTE_GAUNT], gfac;
03657
03658
03659 DEBUG_ENTRY( "PunchGaunts()" );
03660
03661
03662
03663 tempsave = phycon.te;
03664
03665 for( i=0; i < NTE_GAUNT; i++ )
03666 {
03667 ste[i] = 0.5f*i;
03668 }
03669
03670 for( i=0; i < NENR_GAUNT; i++ )
03671 {
03672 ener[i] = 0.5f*i - 8.f;
03673 }
03674
03675 for( charge = 1; charge<LIMELM; charge++ )
03676 {
03677
03678 z = (float)log10((double)charge);
03679
03680
03681 for( ite=0; ite < NTE_GAUNT; ite++ )
03682 {
03683 phycon.alogte = ste[ite];
03684 phycon.te = pow(10.,phycon.alogte);
03685
03686 for( j=0; j < NENR_GAUNT; j++ )
03687 {
03688 gfac = cont_gaunt_calc( phycon.te, (double)charge, pow(10.,(double)ener[j]) );
03689
03690
03691 g[j][ite] = gfac;
03692 }
03693 }
03694
03695
03696 fprintf( ioPUN, " " );
03697 for( i=1; i <= NTE_GAUNT; i++ )
03698 {
03699 fprintf( ioPUN, "\t%6.1f", ste[i-1] );
03700 }
03701 fprintf( ioPUN, "\n" );
03702
03703 for( j=0; j < NENR_GAUNT; j++ )
03704 {
03705 fprintf( ioPUN, "\t%6.2f", ener[j] );
03706 for( ite=0; ite < NTE_GAUNT; ite++ )
03707 {
03708 fprintf( ioPUN, "\t%6.2f", log10(g[j][ite]) );
03709 }
03710 fprintf( ioPUN, "\n" );
03711 }
03712
03713 fprintf( ioPUN, " " );
03714 for( i=0; i < NTE_GAUNT; i++ )
03715 {
03716 fprintf( ioPUN, "\t%6.1f", ste[i] );
03717 }
03718 fprintf( ioPUN, "\n" );
03719
03720
03721 fprintf( ioPUN, " " );
03722 for( i=0; i < NTE_GAUNT; i++ )
03723 {
03724 fprintf( ioPUN, "\t%6.1f", ste[i] );
03725 }
03726 fprintf( ioPUN, "\n" );
03727
03728 for( i=0; i < NTE_GAUNT; i++ )
03729 {
03730 for( j=0; j < NENR_GAUNT; j++ )
03731 {
03732 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", ste[i], ener[j],
03733 log10(g[j][i]) );
03734 }
03735 }
03736
03737 fprintf( ioPUN, "Below is log(u), log(gamma**2), gff\n" );
03738
03739 for( i=0; i < NTE_GAUNT; i++ )
03740 {
03741 for( j=0; j < NENR_GAUNT; j++ )
03742 {
03743 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", 2.*z + log10(TE1RYD) - ste[i] , log10(TE1RYD)+ener[j]-ste[i],
03744 log10(g[j][i]) );
03745 }
03746 }
03747 fprintf( ioPUN, "end of charge = %li\n", charge);
03748 fprintf( ioPUN, "****************************\n");
03749 }
03750
03751 phycon.te = tempsave;
03752 phycon.alogte = log10( phycon.te );
03753
03754 DEBUG_EXIT( "PunchGaunts()" );
03755 return;
03756 }