cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*PunchDo produce punch output during calculation, 00004 * chTime is 'MIDL' during calculation, 'LAST' at the end */ 00005 /*PunchNewContinuum produce the 'punch new continuum' output */ 00006 /*PunchLineStuff punch optical depths or source functions for all transferred lines */ 00007 /*pun1Line called by PunchLineStuff to produce output for one line */ 00008 /*PunchNewContinuum produce the 'punch new continuum' output */ 00009 /*PunLineIntensity produce the 'punch lines intensity' output */ 00010 /* punch h emission, for AGN3 chapter 4, routine is below */ 00011 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */ 00012 /* the number of emission lines across one line of printout */ 00013 /*PunchSpecial generate output for the punch special command */ 00014 /*punResults punch results from punch results command */ 00015 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */ 00016 /*PunchGaunts called by punch gaunts command to output gaunt factors */ 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 "hcmap.h" 00064 #include "input.h" 00065 #include "punch.h" 00066 #include "optimize.h" 00067 #include "grid.h" 00068 00069 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */ 00070 /* the number of emission lines across one line of printout */ 00071 STATIC void PunResults1Line( 00072 FILE * ioPUN, 00073 /* 4 char + null string */ 00074 const char *chLab, 00075 realnum wl, 00076 double xInten, 00077 const char *chFunction);/* null terminated string saying what to do */ 00078 00079 /*PunchGaunts called by punch gaunts command to output gaunt factors */ 00080 STATIC void PunchGaunts(FILE* ioPUN); 00081 00082 /*punResults punch results from punch results command */ 00083 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */ 00084 STATIC void punResults(FILE* ioPUN); 00085 00086 STATIC void PunchLineStuff( 00087 FILE * ioPUN, 00088 const char *chJob , 00089 realnum xLimit); 00090 00091 /* punch h emission, for chapter 4, routine is below */ 00092 STATIC void AGN_Hemis(FILE *ioPUN ); 00093 00094 /*PunchNewContinuum produce the 'punch new continuum' output */ 00095 STATIC void PunchNewContinuum(FILE * ioPUN , long ipCon ); 00096 00097 /*PunLineIntensity produce the 'punch lines intensity' output */ 00098 STATIC void PunLineIntensity(FILE * ioPUN); 00099 00100 char *chDummy; 00101 00102 void PunchDo( 00103 /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */ 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 * the "last" option on punch command, to punch on last iteration, 00147 * is parsed at the top of the loop in only one place. 00148 * no further action is needed at all for punch last to work 00149 * ok throughout this routine 00150 */ 00151 00152 /* 00153 * each branch can have a test whether chTime is or is not "LAST" 00154 * 00155 * if( strcmp(chTime,"LAST") == 0 ) <== print after iteration is complete 00156 * 00157 * if "LAST" then this is last call to routine after iteration complete 00158 * punch only if "LAST" when results at end of iteration are needed 00159 * 00160 * if( strcmp(chTime,"LAST") != 0 ) <== print for every zone 00161 * 00162 * test for .not."LAST" is for every zone result, where you do not 00163 * want to punch last zone twice 00164 */ 00165 00166 /* return if no punch to do */ 00167 if( punch.npunch < 1 ) 00168 { 00169 return; 00170 } 00171 00172 /* during a grid calculation this routine saves grid points after 00173 * cloudy is called. we may output it below */ 00174 if( grid.lgGrid ) 00175 { 00176 if( chTime[0]=='L' ) 00177 GridGatherInCloudy(); 00178 00179 /* save grid information */ 00180 GridGatherAfterCloudy( chTime ); 00181 } 00182 00183 for( ipPun=0; ipPun < punch.npunch; ipPun++ ) 00184 { 00185 /* this global variable to remember where in the punch stack we are */ 00186 punch.ipConPun = ipPun; 00187 00188 /* iterations.lgLastIt is true if this is last iteration 00189 * lgPunLstIter set true if 'last' key occurred on punch command 00190 * normally is false. This will skip punching if last set and 00191 * this is not last iteration */ 00192 if( iterations.lgLastIt || (!punch.lgPunLstIter[ipPun]) ) 00193 { 00194 00195 if( strcmp(punch.chPunch[ipPun],"ABUN") == 0 ) 00196 { 00197 /* punch abundances vs depth */ 00198 if( strcmp(chTime,"LAST") != 0 ) 00199 { 00200 fprintf( punch.ipPnunit[ipPun], "%.2f", 00201 log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) ); 00202 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 00203 { 00204 /* >>chng 05 feb 03, protect against non-positive abundances, 00205 * bug caught by Marcelo Castellanos */ 00206 fprintf( punch.ipPnunit[ipPun], "\t%.2f", 00207 log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) ); 00208 } 00209 fprintf( punch.ipPnunit[ipPun], "\n" ); 00210 } 00211 } 00212 00213 else if( strcmp(punch.chPunch[ipPun],"21CM") == 0 ) 00214 { 00215 /* punch information about 21 cm line */ 00216 if( strcmp(chTime,"LAST") != 0 ) 00217 { 00218 fprintf( punch.ipPnunit[ipPun], 00219 "%.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", 00220 /* depth, cm */ 00221 radius.depth_mid_zone, 00222 hyperfine.Tspin21cm , 00223 phycon.te , 00224 /* temperature from Lya - 21 cm optical depth ratio */ 00225 3.84e-7* Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon / 00226 SDIV( HFLines[0].Emis->TauCon ), 00227 /*TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ),*/ 00228 HFLines[0].Lo->Pop , 00229 HFLines[0].Hi->Pop , 00230 OccupationNumberLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ), 00231 HFLines[0].Emis->TauCon , 00232 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon, 00233 HFLines[0].Emis->PopOpc, 00234 /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */ 00235 (dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop)/SDIV( hyperfine.Tspin21cm), 00236 /* why was above multiplied by this following term? */ 00237 /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/ 00238 HFLines[0].Emis->TauIn, 00239 TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) , 00240 colden.H0_ov_Tspin, 00241 /*>>chng 27 mar, GS, integrated 21cm spin temperature*/ 00242 colden.H0_21cm_lower, 00243 colden.H0_21cm_upper, 00244 -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower) 00245 ); 00246 } 00247 } 00248 00249 else if( strcmp(punch.chPunch[ipPun],"AGES") == 0 ) 00250 { 00251 /* punch timescales vs depth */ 00252 if( strcmp(chTime,"LAST") != 0 ) 00253 { 00254 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00255 /* depth, cm */ 00256 radius.depth_mid_zone, 00257 /* cooling timescale */ 00258 dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot, 00259 /* H2 destruction timescale */ 00260 timesc.time_H2_Dest_here, 00261 /* CO destruction timescale */ 00262 timesc.AgeCOMoleDest[findspecies("CO")->index], 00263 /* OH destruction timescale */ 00264 timesc.AgeCOMoleDest[findspecies("OH")->index], 00265 /* H recombination timescale */ 00266 1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) ); 00267 } 00268 } 00269 00270 else if( strcmp(punch.chPunch[ipPun]," AGN") == 0 ) 00271 { 00272 if( strcmp(chTime,"LAST") == 0 ) 00273 { 00274 if( strcmp( punch.chPunchArgs[ipPun], "HECS" ) == 0 ) 00275 { 00276 /* this routine is in helike.c */ 00277 AGN_He1_CS(punch.ipPnunit[ipPun]); 00278 } 00279 if( strcmp( punch.chPunchArgs[ipPun], "HEMI" ) == 0 ) 00280 { 00281 /* punch h emiss, for chapter 4, routine is below */ 00282 AGN_Hemis(punch.ipPnunit[ipPun]); 00283 } 00284 else 00285 { 00286 fprintf( ioQQQ, " PunchDo does not recognize flag %4.4s set for AGN punch. This is impossible.\n", 00287 punch.chPunch[ipPun] ); 00288 ShowMe(); 00289 cdEXIT(EXIT_FAILURE); 00290 } 00291 } 00292 } 00293 00294 else if( strcmp(punch.chPunch[ipPun],"ASSE") == 0 ) 00295 { 00296 if( strcmp(chTime,"LAST") == 0 ) 00297 { 00298 /* punch the assert output */ 00299 lgCheckAsserts(punch.ipPnunit[ipPun]); 00300 } 00301 } 00302 00303 else if( strcmp(punch.chPunch[ipPun],"AVER") == 0 ) 00304 { 00305 if( strcmp(chTime,"LAST") == 0 ) 00306 { 00307 /* punch the averages output */ 00308 punch_average( ipPun , "PUNCH", chDummy ); 00309 } 00310 } 00311 00312 else if( strncmp(punch.chPunch[ipPun],"CHA",3) == 0 ) 00313 { 00314 if( strcmp(chTime,"LAST") == 0 ) 00315 { 00316 /* one of the charge transfer options, all in chargtran.c */ 00317 ChargTranPun( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] ); 00318 } 00319 } 00320 00321 else if( strcmp(punch.chPunch[ipPun],"COOL") == 0 ) 00322 { 00323 /* punch cooling, routine in file of same name */ 00324 if( strcmp(chTime,"LAST") != 0 ) 00325 CoolPunch(punch.ipPnunit[ipPun]); 00326 } 00327 00328 else if( strncmp(punch.chPunch[ipPun],"DYN" , 3) == 0 ) 00329 { 00330 /* punch dynamics xxx, information about dynamical solutions */ 00331 if( strcmp(chTime,"LAST") != 0 ) 00332 DynaPunch(punch.ipPnunit[ipPun] ,punch.chPunch[ipPun][3] ); 00333 } 00334 00335 else if( strcmp(punch.chPunch[ipPun],"ENTH") == 0 ) 00336 { 00337 if( strcmp(chTime,"LAST") != 0 ) 00338 fprintf( punch.ipPnunit[ipPun], 00339 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00340 radius.depth_mid_zone, 00341 phycon.EnthalpyDensity, 00342 phycon.EnergyExcitation, 00343 phycon.EnergyIonization, 00344 phycon.EnergyBinding , 00345 0.5*POW2(wind.windv)*dense.xMassDensity , /* KE */ 00346 5./2.*pressure.PresGasCurr , /* thermal plus PdV work */ 00347 magnetic.EnthalpyDensity); /* magnetic terms */ 00348 } 00349 00350 else if( strcmp(punch.chPunch[ipPun],"COLU") == 0 ) 00351 { 00352 /* punch column densities */ 00353 if( strcmp(chTime,"LAST") == 0 ) 00354 { 00355 PrtColumns(punch.ipPnunit[ipPun]); 00356 } 00357 } 00358 00359 else if( strcmp(punch.chPunch[ipPun],"COLS") == 0 ) 00360 { 00361 /* punch some column densities */ 00362 if( strcmp(chTime,"LAST") == 0 ) 00363 { 00364 char chHeader[2]; 00365 punch_colden(chHeader , punch.ipPnunit[ipPun] , "PUNS" ); 00366 } 00367 } 00368 00369 else if( strcmp(punch.chPunch[ipPun],"COMP") == 0 ) 00370 { 00371 /* Compton energy exchange coefficients */ 00372 if( strcmp(chTime,"LAST") != 0 ) 00373 { 00374 for( jj=0; jj<rfield.nflux; jj = jj + punch.ncPunchSkip) 00375 { 00376 fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e\n", 00377 rfield.anu[jj], rfield.comup[jj]/rfield.widflx[jj], 00378 rfield.comdn[jj]/rfield.widflx[jj] ); 00379 } 00380 } 00381 } 00382 00383 /* punch continuum commands */ 00384 else if( strcmp(punch.chPunch[ipPun],"CON ") == 0 ) 00385 { 00386 /* this is the usual "punch continuum" case */ 00387 /* >>chng 06 apr 03, add every option to do every zone */ 00388 /* if punch every is set then nPunchEveryZone must be positive 00389 * was init to -1 */ 00390 bool lgPrintThis =false; 00391 if( punch.lgPunchEveryZone[ipPun] ) 00392 { 00393 /* this branch, every option is on line so want to print every n zone */ 00394 if( strcmp(chTime,"LAST") != 0 ) 00395 { 00396 /* not last zone - print first and intermediate cases */ 00397 if( nzone==1 ) 00398 { 00399 lgPrintThis = true; 00400 } 00401 else if( nzone%punch.nPunchEveryZone[ipPun]==0 ) 00402 { 00403 lgPrintThis = true; 00404 } 00405 } 00406 else 00407 { 00408 /* this is last zone, print only if did not trip on above */ 00409 if( nzone!=1 && nzone%punch.nPunchEveryZone[ipPun]!=0 ) 00410 { 00411 lgPrintThis = true; 00412 } 00413 } 00414 } 00415 else 00416 { 00417 /* this branch, not "every", so only print the last zone */ 00418 if( strcmp(chTime,"LAST") == 0 ) 00419 lgPrintThis = true; 00420 } 00421 ASSERT( !punch.lgPunchEveryZone[ipPun] || punch.nPunchEveryZone[ipPun]>0 ); 00422 if( lgPrintThis ) 00423 { 00424 if( punch.lgPunchEveryZone[ipPun] && nzone!=1) 00425 fprintf( punch.ipPnunit[ipPun], "%s\n", 00426 punch.chHashString ); 00427 00428 for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip) 00429 { 00430 /* four continua predicted here; 00431 * incident, attenuated incident, emitted, 00432 * then attenuated incident + emitted, last reflected 00433 * reflected continuum is stored relative to inner radius 00434 * others are stored for this radius */ 00435 00436 /* NB this code also used in punch emitted, 00437 * transmitted continuum commands */ 00438 00439 /* the incident continuum */ 00440 flxin = rfield.flux_total_incident[j]*rfield.anu2[j]* 00441 EN1RYD/rfield.widflx[j]; 00442 00443 /* the reflected continuum */ 00444 flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/rfield.widflx[j] + 00445 rfield.anu[j]*punch.PunchLWidth*rfield.reflin[j])*EN1RYD; 00446 00447 /* the attenuated incident continuum */ 00448 flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/ 00449 rfield.widflx[j]*radius.r1r0sq * rfield.trans_coef_total[j]; 00450 00451 /* the outward emitted continuum */ 00452 conem = ((rfield.ConEmitOut[j])/ 00453 rfield.widflx[j]*rfield.anu2[j] + punch.PunchLWidth* 00454 rfield.outlin[j]*rfield.anu[j])*radius.r1r0sq* 00455 EN1RYD*geometry.covgeo; 00456 00457 /* sum of emitted and transmitted continua */ 00458 flxtrn = conem + flxatt; 00459 00460 /* photon energy in appropriate energy or wavelength units */ 00461 fprintf( punch.ipPnunit[ipPun],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) ); 00462 /* incident continuum */ 00463 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxin ); 00464 /* trans cont */ 00465 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxatt ); 00466 /* DiffOut cont */ 00467 fprintf( punch.ipPnunit[ipPun],"%.3e\t", conem ); 00468 /* net trans cont */ 00469 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn ); 00470 /* reflected cont */ 00471 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref ); 00472 /* total cont */ 00473 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref + flxtrn ); 00474 fprintf( punch.ipPnunit[ipPun], "%4.4s\t%4.4s\t", 00475 /* line label */ 00476 rfield.chLineLabel[j] , 00477 /* cont label*/ 00478 rfield.chContLabel[j] ); 00479 /* number of lines within that cell over cell width 00480 * punch raw continuum has number by itself */ 00481 fprintf( punch.ipPnunit[ipPun], "%.2f\n", rfield.line_count[j]/rfield.widflx[j]*rfield.anu[j] ); 00482 } 00483 } 00484 } 00485 00486 /* this is the punch spectrum command, 00487 * the new continuum command that will replace the previous one */ 00488 else if( strcmp(punch.chPunch[ipPun],"CONN") == 0 ) 00489 { 00490 if( strcmp(chTime,"LAST") == 0 ) 00491 /* io unit and which new continuum this is (was set when punch read in */ 00492 PunchNewContinuum( punch.ipPnunit[ipPun] , (long)punch.punarg[ipPun][0]); 00493 } 00494 00495 else if( strcmp(punch.chPunch[ipPun],"CONC") == 0 ) 00496 { 00497 /* punch incident continuum */ 00498 /* set pointer for possible change in units of energy in continuum 00499 * AnuUnit will give anu in whatever units were set with punch units */ 00500 if( strcmp(chTime,"LAST") == 0 ) 00501 { 00502 /* incident continuum */ 00503 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip) 00504 { 00505 flxin = rfield.flux_total_incident[j]*rfield.anu2[j]* 00506 EN1RYD/rfield.widflx[j]; 00507 /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */ 00508 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\n", 00509 AnuUnit(rfield.AnuOrg[j]), flxin ); 00510 } 00511 } 00512 } 00513 00514 else if( strcmp(punch.chPunch[ipPun],"CONG") == 0 ) 00515 { 00516 /* punch emitted grain continuum in optically thin limit */ 00517 if( strcmp(chTime,"LAST") == 0 ) 00518 { 00519 /* GrainMakeDiffuse broke out emission into types 00520 * according to matType */ 00521 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip) 00522 { 00523 fsum = gv.GraphiteEmission[j]*rfield.anu2[j]* 00524 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo; 00525 00526 fout = gv.SilicateEmission[j]*rfield.anu2[j]* 00527 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo; 00528 00529 /* anu is .5e format to resolve energy mesh near 1 Ryd 00530 * AnuUnit gives anu in whatever units were set with units option */ 00531 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n", 00532 AnuUnit(rfield.AnuOrg[j]) , fsum , fout , 00533 /* total emission per unit volume defined in GrainMakeDiffuse 00534 * used in RT_diffuse to form total grain emission */ 00535 gv.GrainEmission[j]*rfield.anu2[j]* 00536 EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo ); 00537 } 00538 } 00539 } 00540 00541 else if( strcmp(punch.chPunch[ipPun],"CONR") == 0 ) 00542 { 00543 /* punch reflected continuum */ 00544 /* set pointer for possible change in units of energy in continuum 00545 * AnuUnit will give anu in whatever units were set with punch units */ 00546 if( strcmp(chTime,"LAST") == 0 ) 00547 { 00548 if( geometry.lgSphere ) 00549 { 00550 fprintf( punch.ipPnunit[ipPun], " Reflected continuum not predicted when SPHERE is set.\n" ); 00551 fprintf( ioQQQ , 00552 "\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" ); 00553 cdEXIT(EXIT_FAILURE); 00554 } 00555 00556 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip) 00557 { 00558 /* reflected continuum */ 00559 flxref = rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/ 00560 rfield.widflx[j]*EN1RYD; 00561 /* reflected lines */ 00562 fref = rfield.anu[j]*punch.PunchLWidth* 00563 rfield.reflin[j]*EN1RYD; 00564 /* ratio of reflected to incident continuum, the albedo */ 00565 if( rfield.flux_total_incident[j] > 1e-25 ) 00566 { 00567 av = rfield.ConRefIncid[j]/rfield.flux_total_incident[j]; 00568 } 00569 else 00570 { 00571 av = 0.; 00572 } 00573 /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */ 00574 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4s\n", 00575 AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref, 00576 av, rfield.chContLabel[j] ); 00577 } 00578 } 00579 } 00580 00581 else if( strcmp(punch.chPunch[ipPun],"CNVE") == 0 ) 00582 { 00583 /* the punch convergence error command */ 00584 if( strcmp(chTime,"LAST") != 0 ) 00585 { 00586 fprintf( punch.ipPnunit[ipPun], 00587 "%.5e\t%li\t%.4e\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n", 00588 radius.depth_mid_zone, 00589 conv.nPres2Ioniz, 00590 pressure.PresTotlCorrect, 00591 pressure.PresTotlCurr, 00592 (pressure.PresTotlCorrect - pressure.PresTotlCurr)*100./pressure.PresTotlCorrect, 00593 dense.EdenTrue, 00594 dense.eden, 00595 (dense.EdenTrue - dense.eden)*100./dense.EdenTrue, 00596 thermal.htot, 00597 thermal.ctot, 00598 (thermal.htot - thermal.ctot)*100./thermal.htot ); 00599 } 00600 } 00601 00602 else if( strcmp(punch.chPunch[ipPun],"CONB") == 0 ) 00603 { 00604 /* punch continuum bins binning */ 00605 /* set pointer for possible change in units of energy in continuum 00606 * AnuUnit will give anu in whatever units were set with punch units */ 00607 if( strcmp(chTime,"LAST") != 0 ) 00608 { 00609 for( j=0; j<rfield.nupper; j = j + punch.ncPunchSkip) 00610 { 00611 fprintf( punch.ipPnunit[ipPun], "%14.5e %14.5e %14.5e\n", 00612 AnuUnit(rfield.AnuOrg[j]), rfield.anu[j], rfield.widflx[j] ); 00613 } 00614 } 00615 } 00616 00617 else if( strcmp(punch.chPunch[ipPun],"COND") == 0 ) 00618 { 00619 /* punch diffuse continuum the local thermal emission */ 00620 /* set pointer for possible change in units of energy in continuum 00621 * AnuUnit will give anu in whatever units were set with punch units */ 00622 if( strcmp(chTime,"LAST") == 0 ) 00623 { 00624 /* this option to punch diffuse continuum */ 00625 for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip) 00626 { 00627 /* >>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */ 00628 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\n", 00629 AnuUnit(rfield.AnuOrg[j]), 00630 rfield.ConEmitLocal[nzone][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]); 00631 } 00632 } 00633 } 00634 00635 else if( strcmp(punch.chPunch[ipPun],"CONE") == 0 ) 00636 { 00637 /* punch emitted continuum */ 00638 /* set pointer for possible change in units of energy in continuum 00639 * AnuUnit will give anu in whatever units were set with punch units */ 00640 if( strcmp(chTime,"LAST") == 0 ) 00641 { 00642 /* punch emitted continuum */ 00643 for( j=0; j<rfield.nflux;j = j +punch.ncPunchSkip) 00644 { 00645 /* this is the reflected component */ 00646 flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/ 00647 rfield.widflx[j] + rfield.anu[j]*punch.PunchLWidth* 00648 rfield.reflin[j])*EN1RYD; 00649 00650 /* this is the total emission in the outward direction */ 00651 conem = (rfield.ConEmitOut[j])/rfield.widflx[j]*rfield.anu2[j] + 00652 punch.PunchLWidth*rfield.outlin[j]*rfield.anu[j]; 00653 00654 conem *= radius.r1r0sq*EN1RYD*geometry.covgeo; 00655 00656 /* output: photon energy, reflected, outward, total emission 00657 * >>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */ 00658 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n", 00659 AnuUnit(rfield.AnuOrg[j]), 00660 flxref, 00661 conem, 00662 flxref + conem, 00663 rfield.chLineLabel[j], 00664 rfield.chContLabel[j] 00665 ); 00666 } 00667 } 00668 } 00669 00670 /* punch fine continuum command */ 00671 else if( strcmp(punch.chPunch[ipPun],"CONf") == 0 ) 00672 { 00673 if( strcmp(chTime,"LAST") == 0 ) 00674 { 00675 long nu_hi , nskip; 00676 if( punch.punarg[ipPun][0] > 0. ) 00677 /* possible lower bounds to energy range - 00678 * 0 if not set with range option*/ 00679 j = ipFineCont( punch.punarg[ipPun][0] ); 00680 else 00681 j = 0; 00682 00683 /* upper limit set with range option */ 00684 if( punch.punarg[ipPun][1]> 0. ) 00685 nu_hi = ipFineCont( punch.punarg[ipPun][1]); 00686 else 00687 nu_hi = rfield.nfine; 00688 00689 /* number of cells to bring together, default is 10 */ 00690 nskip = (long)punch.punarg[ipPun][2]; 00691 00692 do 00693 { 00694 realnum sum1 = rfield.fine_opt_depth[j]; 00695 realnum xnu = rfield.fine_anu[j]; 00696 for( jj=1; jj<nskip; ++jj ) 00697 { 00698 xnu += rfield.fine_anu[j+jj]; 00699 sum1 += rfield.fine_opt_depth[j+jj]; 00700 } 00701 fprintf( punch.ipPnunit[ipPun], 00702 "%.6e\t%.3e\n", 00703 AnuUnit(xnu/nskip), 00704 sexp(sum1/nskip) ); 00705 j += nskip; 00706 } while( j < nu_hi ); 00707 } 00708 } 00709 00710 else if( strcmp(punch.chPunch[ipPun],"CONi") == 0 ) 00711 { 00712 /* punch continuum interactions */ 00713 /* set pointer for possible change in units of energy in continuum 00714 * AnuUnit will give anu in whatever units were set with punch units */ 00715 00716 /* continuum interactions */ 00717 if( strcmp(chTime,"LAST") != 0 ) 00718 { 00719 /* this is option to set lowest energy */ 00720 if( punch.punarg[ipPun][0] <= 0. ) 00721 { 00722 i1 = 1; 00723 } 00724 else if( punch.punarg[ipPun][0] < 100. ) 00725 { 00726 i1 = ipoint(punch.punarg[ipPun][0]); 00727 } 00728 else 00729 { 00730 i1 = (long int)punch.punarg[ipPun][0]; 00731 } 00732 00733 fref = 0.; 00734 fout = 0.; 00735 fsum = 0.; 00736 sum = 0.; 00737 flxin = 0.; 00738 00739 for( j=i1-1; j < rfield.nflux; j++ ) 00740 { 00741 fref += rfield.flux[j]*opac.opacity_abs[j]; 00742 fout += rfield.otslin[j]*opac.opacity_abs[j]; 00743 fsum += rfield.otscon[j]*opac.opacity_abs[j]; 00744 sum += rfield.ConInterOut[j]*opac.opacity_abs[j]; 00745 flxin += (rfield.outlin[j] + rfield.outlin_noplot[j])*opac.opacity_abs[j]; 00746 } 00747 fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e%10.2e%10.2e\n", 00748 fref, fout, fsum, sum, flxin ); 00749 } 00750 } 00751 00752 else if( strcmp(punch.chPunch[ipPun],"CONI") == 0 ) 00753 { 00754 /* punch ionizing continuum */ 00755 /* set pointer for possible change in units of energy in continuum 00756 * AnuUnit will give anu in whatever units were set with punch units */ 00757 00758 if( (punch.punarg[ipPun][2]>0) || (strcmp(chTime,"LAST") == 0) ) 00759 { 00760 /* this flag will remember whether we have ever printed anything */ 00761 bool lgPrt=false; 00762 if( punch.punarg[ipPun][2]>0 ) 00763 fprintf(punch.ipPnunit[ipPun],"#punch every zone %li\n", nzone); 00764 00765 /* punch ionizing continuum command 00766 * this is option to set lowest energy, 00767 * if no number was entered then this was zero */ 00768 if( punch.punarg[ipPun][0] <= 0. ) 00769 { 00770 i1 = 1; 00771 } 00772 else if( punch.punarg[ipPun][0] < 100. ) 00773 { 00774 i1 = ipoint(punch.punarg[ipPun][0]); 00775 } 00776 else 00777 { 00778 i1 = (long int)punch.punarg[ipPun][0]; 00779 } 00780 00781 sum = 0.; 00782 for( j=i1-1; j < rfield.nflux; j++ ) 00783 { 00784 flxcor = rfield.flux[j] + 00785 rfield.otslin[j] + 00786 rfield.otscon[j] + 00787 rfield.ConInterOut[j] + 00788 rfield.outlin[j] + rfield.outlin_noplot[j]; 00789 00790 sum += flxcor*opac.opacity_abs[j]; 00791 } 00792 00793 if( sum > 0. ) 00794 { 00795 sum = 1./sum; 00796 } 00797 else 00798 { 00799 sum = 1.; 00800 } 00801 00802 fsum = 0.; 00803 00804 for( j=i1-1; j<rfield.nflux; ++j) 00805 { 00806 flxcor = rfield.flux[j] + 00807 rfield.otslin[j] + 00808 rfield.otscon[j] + 00809 rfield.ConInterOut[j]+ 00810 rfield.outlin[j] + rfield.outlin_noplot[j]; 00811 00812 fsum += flxcor*opac.opacity_abs[j]; 00813 00814 /* punched quantities are freq, flux, flux*cross sec, 00815 * fraction of total, integral fraction of total */ 00816 RateInter = flxcor*opac.opacity_abs[j]*sum; 00817 00818 /* punage(ipPun,2) is lowest interaction rate to consider, def=0.01 (1 percent) */ 00819 /* >>chng 01 nov 22, format to c-friendly */ 00820 if( (RateInter >= punch.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) ) 00821 { 00822 lgPrt = true; 00823 /* >>chng 96 oct 22, format of anu to 11.5 to resolve energy mesh near 1 Ryd */ 00824 fprintf( punch.ipPnunit[ipPun], 00825 "%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", 00826 j, 00827 AnuUnit(rfield.AnuOrg[j]), 00828 flxcor, 00829 flxcor*opac.opacity_abs[j], 00830 rfield.flux[j]/flxcor, 00831 rfield.otslin[j]/flxcor, 00832 rfield.otscon[j]/flxcor, 00833 (rfield.outlin[j] + rfield.outlin_noplot[j])/flxcor, 00834 rfield.ConInterOut[j]/flxcor, 00835 RateInter, 00836 fsum*sum, 00837 rfield.chLineLabel[j], 00838 rfield.chContLabel[j] ); 00839 } 00840 } 00841 if( !lgPrt ) 00842 { 00843 /* entered logical block but did not print anything */ 00844 fprintf(punch.ipPnunit[ipPun], 00845 " punchdo, the PUNCH IONIZING CONTINUUM command " 00846 "did not find a strong point, sum and fsum were %.2e %.2e\n", 00847 sum,fsum); 00848 fprintf(punch.ipPnunit[ipPun], 00849 " punchdo, the low-frequency energy was %.5e Ryd\n", 00850 rfield.anu[i1-1]); 00851 fprintf(punch.ipPnunit[ipPun], 00852 " You can reset the threshold for the lowest fractional " 00853 "interaction to print with the second number of the punch command\n" 00854 " The fraction was %.3f and this was too large.\n", 00855 punch.punarg[ipPun][1]); 00856 } 00857 } 00858 } 00859 00860 else if( strcmp(punch.chPunch[ipPun],"CORA") == 0 ) 00861 { 00862 /* punch raw continuum */ 00863 /* set pointer for possible change in units of energy in continuum 00864 * AnuUnit will give anu in whatever units were set with punch units */ 00865 00866 if( strcmp(chTime,"LAST") == 0 ) 00867 { 00868 /* this option to punch all raw ionizing continuum */ 00869 for( j=0;j<rfield.nflux;j = j + punch.ncPunchSkip) 00870 { 00871 fprintf( punch.ipPnunit[ipPun], 00872 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t", 00873 AnuUnit(rfield.AnuOrg[j]), 00874 rfield.flux[j], 00875 rfield.otslin[j], 00876 rfield.otscon[j], 00877 rfield.ConRefIncid[j], 00878 rfield.ConEmitReflec[j], 00879 rfield.ConInterOut[j], 00880 rfield.outlin[j]+rfield.outlin_noplot[j], 00881 rfield.ConEmitOut[j] , 00882 rfield.chLineLabel[j], 00883 rfield.chContLabel[j] 00884 ); 00885 /* number of lines within that cell */ 00886 fprintf( punch.ipPnunit[ipPun], "%li\n", rfield.line_count[j] ); 00887 } 00888 } 00889 } 00890 00891 else if( strcmp(punch.chPunch[ipPun],"CONo") == 0 ) 00892 { 00893 /* punch outward local continuum */ 00894 /* set pointer for possible change in units of energy in continuum 00895 * AnuUnit will give anu in whatever units were set with punch units */ 00896 00897 if( strcmp(chTime,"LAST") == 0 ) 00898 { 00899 /* option to punch out outward continuum here */ 00900 for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip) 00901 { 00902 fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e\n", 00903 AnuUnit(rfield.AnuOrg[j]), 00904 rfield.ConEmitOut[j]+ rfield.outlin[j] + rfield.outlin_noplot[j], 00905 (rfield.flux[j] + rfield.otscon[j] + rfield.otslin[j] + 00906 rfield.ConInterOut[j])*opac.opacity_abs[j]* 00907 rfield.anu[j] ); 00908 } 00909 } 00910 } 00911 00912 else if( strcmp(punch.chPunch[ipPun],"CONO") == 0 ) 00913 { 00914 /* punch outward continuum */ 00915 /* set pointer for possible change in units of energy in continuum 00916 * AnuUnit will give anu in whatever units were set with punch units */ 00917 00918 if( strcmp(chTime,"LAST") == 0 ) 00919 { 00920 /* option to punch out outward continuum */ 00921 for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip) 00922 { 00923 fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e%10.2e%10.2e\n", 00924 AnuUnit(rfield.AnuOrg[j]), 00925 rfield.flux[j]*rfield.anu2[j]* EN1RYD/rfield.widflx[j]*radius.r1r0sq, 00926 rfield.ConInterOut[j]/rfield.widflx[j]*rfield.anu2[j]*radius.r1r0sq*EN1RYD, 00927 punch.PunchLWidth* (rfield.outlin[j]+rfield.outlin_noplot[j])*rfield.anu[j]*radius.r1r0sq*EN1RYD, 00928 (rfield.ConInterOut[j]/rfield.widflx[j]* 00929 rfield.anu2[j] + punch.PunchLWidth*(rfield.outlin[j]+rfield.outlin_noplot[j])* 00930 rfield.anu[j])*radius.r1r0sq*EN1RYD ); 00931 } 00932 } 00933 } 00934 00935 else if( strcmp(punch.chPunch[ipPun],"CONT") == 0 ) 00936 { 00937 /* punch transmitted continuum - this is not the main "punch continuum" 00938 * command - search on "CON " above 00939 * set pointer for possible change in units of energy in continuum 00940 * AnuUnit will give anu in whatever units were set with punch units */ 00941 00942 if( strcmp(chTime,"LAST") == 0 ) 00943 { 00944 /* this option to punch transmitted continuum */ 00945 for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip) 00946 { 00947 /* attenuated incident continuum 00948 * >>chng 97 jul 10, remove PunchLWidth from this one only since 00949 * we must conserve energy even in lines 00950 * >>chng 07 apr 26 include transmission coefficient */ 00951 flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/ 00952 rfield.widflx[j]*radius.r1r0sq * rfield.trans_coef_total[j]; 00953 00954 /*conem = (rfield.ConOutNoInter[j] + rfield.ConInterOut[j]+rfield.outlin[j])* 00955 rfield.anu2[j]; 00956 conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;*/ 00957 /* >>chng 00 jan 03, above did not include all contributors. 00958 * Pasted in below from usual 00959 * punch continuum command */ 00960 /* >>chng 04 jul 15, removed factor of punch.PunchLWidth - 00961 * this should not be there to conserve energy, as explained in hazy 00962 * where command was documented, and in comment above. caught by PvH */ 00963 /* >>chng 04 jul 23, incorrect use of outlin - before multiplied by an2, 00964 * quantity should be photons per Ryd, since init quantity is 00965 * photons per cell. Must div by widflx. caught by PvH */ 00966 /*conem = (rfield.ConEmitOut[j]/rfield.widflx[j]*rfield.anu2[j] + 00967 rfield.outlin[j]*rfield.anu[j])*radius.r1r0sq* 00968 EN1RYD*geometry.covgeo;*/ 00969 conem = (rfield.ConEmitOut[j] + rfield.outlin[j]) / rfield.widflx[j]* 00970 rfield.anu2[j]*radius.r1r0sq*EN1RYD*geometry.covgeo; 00971 00972 flxtrn = conem + flxatt; 00973 00974 /* use AnuOrg here instead of anu since probably 00975 * going to be used by table read 00976 * and we want good anu array for sanity check*/ 00977 fprintf( punch.ipPnunit[ipPun],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) ); 00978 fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn); 00979 fprintf( punch.ipPnunit[ipPun],"%.3e\n", rfield.trans_coef_total[j] ); 00980 } 00981 } 00982 } 00983 00984 else if( strcmp(punch.chPunch[ipPun],"CON2") == 0 ) 00985 { 00986 /* punch total two-pohoton continuum */ 00987 if( strcmp(chTime,"LAST") == 0 ) 00988 { 00989 /* this option to punch diffuse continuum */ 00990 for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip) 00991 { 00992 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\n", 00993 AnuUnit(rfield.AnuOrg[j]), 00994 rfield.TotDiff2Pht[j]/rfield.widflx[j] , 00995 rfield.TotDiff2Pht[j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]); 00996 } 00997 } 00998 } 00999 01000 else if( strcmp(punch.chPunch[ipPun],"DUSE") == 0 ) 01001 { 01002 /* punch grain extinction - includes only grain opacity, not total */ 01003 if( strcmp(chTime,"LAST") != 0 ) 01004 { 01005 fprintf( punch.ipPnunit[ipPun], " %.5e\t", 01006 radius.depth_mid_zone ); 01007 01008 /* visual extinction of an extended source (like a PDR)*/ 01009 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended); 01010 01011 /* visual extinction of point source (star)*/ 01012 fprintf( punch.ipPnunit[ipPun], "%.2e\n" , rfield.extin_mag_V_point); 01013 } 01014 } 01015 01016 else if( strcmp(punch.chPunch[ipPun],"DUSO") == 0 ) 01017 { 01018 /* punch grain opacity */ 01019 if( strcmp(chTime,"LAST") == 0 ) 01020 { 01021 for( j=0; j < rfield.nflux; j++ ) 01022 { 01023 double scat; 01024 fprintf( punch.ipPnunit[ipPun], 01025 "%.5e\t%.2e\t%.2e\t%.2e\t", 01026 /* photon energy or wavelength */ 01027 AnuUnit(rfield.AnuOrg[j]), 01028 /* total opacity, discount forward scattering */ 01029 gv.dstab[j] + gv.dstsc[j], 01030 /* absorption opacity */ 01031 gv.dstab[j], 01032 /* scatter, with forward discounted */ 01033 gv.dstsc[j] ); 01034 /* add together total scattering, discounting 1-g */ 01035 scat = 0.; 01036 /* sum over all grain species */ 01037 for( nd=0; nd < gv.nBin; nd++ ) 01038 { 01039 scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund; 01040 } 01041 /* finally, scattering including effects of forward scattering */ 01042 fprintf( punch.ipPnunit[ipPun], 01043 "%.2e", scat ); 01044 fprintf( punch.ipPnunit[ipPun], 01045 "%.2e\n", gv.dstsc[j]/(gv.dstab[j] + gv.dstsc[j]) ); 01046 } 01047 } 01048 } 01049 01050 /* punch grain abundance and punch grain D/G ratio commands */ 01051 else if( strcmp(punch.chPunch[ipPun],"DUSA") == 0 || 01052 strcmp(punch.chPunch[ipPun],"DUSD") == 0 ) 01053 { 01054 bool lgDGRatio = ( strcmp(punch.chPunch[ipPun],"DUSD") == 0 ); 01055 01056 /* grain abundance */ 01057 if( strcmp(chTime,"LAST") != 0 ) 01058 { 01059 /* used to print header exactly one time */ 01060 static bool lgMustPrtHeaderDRRatio = true, 01061 lgMustPrtHeaderAbundance=true; 01062 /* print grain headder first if this has not yet been done */ 01063 if( ( lgMustPrtHeaderDRRatio && lgDGRatio ) || 01064 ( lgMustPrtHeaderAbundance && !lgDGRatio ) ) 01065 { 01066 /* only print one header for each case, but must 01067 * track separately if both used in same sim */ 01068 if( lgMustPrtHeaderDRRatio && lgDGRatio ) 01069 lgMustPrtHeaderDRRatio = false; 01070 else if( lgMustPrtHeaderAbundance &&!lgDGRatio ) 01071 lgMustPrtHeaderAbundance = false; 01072 01073 fprintf( punch.ipPnunit[ipPun], "#Depth" ); 01074 for( nd=0; nd < gv.nBin; ++nd ) 01075 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); 01076 fprintf( punch.ipPnunit[ipPun], "\ttotal\n" ); 01077 01078 /* now print grain radius converting from cm to microns */ 01079 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" ); 01080 for( nd=0; nd < gv.nBin; ++nd ) 01081 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); 01082 fprintf( punch.ipPnunit[ipPun], "\txxxx\n" ); 01083 } 01084 fprintf( punch.ipPnunit[ipPun], " %.5e", 01085 radius.depth_mid_zone ); 01086 /* grain abundance per bin in g/cm^3 */ 01087 double total = 0.; 01088 for( nd=0; nd < gv.nBin; ++nd ) 01089 { 01090 double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]* 01091 gv.bin[nd]->cnv_H_pCM3; 01092 if( lgDGRatio ) 01093 abund /= dense.xMassDensity; 01094 fprintf( punch.ipPnunit[ipPun], "\t%.3e", abund ); 01095 total += abund; 01096 } 01097 fprintf( punch.ipPnunit[ipPun], "\t%.3e\n", total ); 01098 } 01099 } 01100 01101 else if( strcmp(punch.chPunch[ipPun],"DUSP") == 0 ) 01102 { 01103 /* grain potential */ 01104 if( strcmp(chTime,"LAST") != 0 ) 01105 { 01106 /* used to print header exactly one time */ 01107 static bool lgMustPrtHeader = true; 01108 /* do labels first if this is first zone */ 01109 if( lgMustPrtHeader ) 01110 { 01111 /* first print string giving grain id */ 01112 fprintf( punch.ipPnunit[ipPun], "#Depth" ); 01113 for( nd=0; nd < gv.nBin; ++nd ) 01114 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); 01115 fprintf( punch.ipPnunit[ipPun], "\n" ); 01116 01117 /* now print grain radius converting from cm to microns */ 01118 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" ); 01119 for( nd=0; nd < gv.nBin; ++nd ) 01120 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); 01121 fprintf( punch.ipPnunit[ipPun], "\n" ); 01122 01123 /* don't need to do this, ever again */ 01124 lgMustPrtHeader = false; 01125 } 01126 fprintf( punch.ipPnunit[ipPun], " %.5e", 01127 radius.depth_mid_zone ); 01128 /* grain potential in eV */ 01129 for( nd=0; nd < gv.nBin; ++nd ) 01130 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD ); 01131 fprintf( punch.ipPnunit[ipPun], "\n" ); 01132 } 01133 } 01134 01135 else if( strcmp(punch.chPunch[ipPun],"DUSR") == 0 ) 01136 { 01137 /* grain H2 formation rates */ 01138 if( strcmp(chTime,"LAST") != 0 ) 01139 { 01140 /* used to print header exactly one time */ 01141 static bool lgMustPrtHeader = true; 01142 /* do labels first if this is first zone */ 01143 if( lgMustPrtHeader ) 01144 { 01145 /* first print string giving grain id */ 01146 fprintf( punch.ipPnunit[ipPun], "#Depth" ); 01147 for( nd=0; nd < gv.nBin; ++nd ) 01148 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); 01149 fprintf( punch.ipPnunit[ipPun], "\n" ); 01150 01151 /* now print grain radius converting from cm to microns */ 01152 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" ); 01153 for( nd=0; nd < gv.nBin; ++nd ) 01154 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); 01155 fprintf( punch.ipPnunit[ipPun], "\n" ); 01156 01157 /* don't need to do this, ever again */ 01158 lgMustPrtHeader = false; 01159 } 01160 fprintf( punch.ipPnunit[ipPun], " %.5e", 01161 radius.depth_mid_zone ); 01162 /* grain formation rate for H2 */ 01163 for( nd=0; nd < gv.nBin; ++nd ) 01164 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used ); 01165 fprintf( punch.ipPnunit[ipPun], "\n" ); 01166 } 01167 } 01168 01169 else if( strcmp(punch.chPunch[ipPun],"DUST") == 0 ) 01170 { 01171 /* grain temperatures - K*/ 01172 if( strcmp(chTime,"LAST") != 0 ) 01173 { 01174 /* used to print header exactly one time */ 01175 static bool lgMustPrtHeader = true; 01176 /* do labels first if this is first zone */ 01177 if( lgMustPrtHeader ) 01178 { 01179 /* first print string giving grain id */ 01180 fprintf( punch.ipPnunit[ipPun], "#Depth" ); 01181 for( nd=0; nd < gv.nBin; ++nd ) 01182 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); 01183 fprintf( punch.ipPnunit[ipPun], "\n" ); 01184 01185 /* now print grain radius converting from cm to microns */ 01186 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" ); 01187 for( nd=0; nd < gv.nBin; ++nd ) 01188 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); 01189 fprintf( punch.ipPnunit[ipPun], "\n" ); 01190 01191 /* don't need to do this, ever again */ 01192 lgMustPrtHeader = false; 01193 } 01194 fprintf( punch.ipPnunit[ipPun], " %.5e", 01195 radius.depth_mid_zone ); 01196 for( nd=0; nd < gv.nBin; ++nd ) 01197 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust ); 01198 fprintf( punch.ipPnunit[ipPun], "\n" ); 01199 } 01200 } 01201 01202 else if( strcmp(punch.chPunch[ipPun],"DUSC") == 0 ) 01203 { 01204 /* punch grain charge - eden from grains and 01205 * charge per grain in electrons / grain */ 01206 if( strcmp(chTime,"LAST") != 0 ) 01207 { 01208 /* used to print header exactly one time */ 01209 static bool lgMustPrtHeader = true; 01210 /* do labels first if this is first zone */ 01211 if( lgMustPrtHeader ) 01212 { 01213 /* first print string giving grain id */ 01214 fprintf( punch.ipPnunit[ipPun], "#Depth\tne(grn)" ); 01215 for( nd=0; nd < gv.nBin; ++nd ) 01216 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); 01217 fprintf( punch.ipPnunit[ipPun], "\n" ); 01218 01219 /* now print grain radius converting from cm to microns */ 01220 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)\tne\t" ); 01221 for( nd=0; nd < gv.nBin; ++nd ) 01222 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); 01223 fprintf( punch.ipPnunit[ipPun], "\n" ); 01224 01225 /* don't need to do this, ever again */ 01226 lgMustPrtHeader = false; 01227 } 01228 01229 fprintf( punch.ipPnunit[ipPun], " %.5e\t%.4e", 01230 radius.depth_mid_zone , 01231 /* electron density contributed by grains, in e/cm^3, 01232 * positive number means grain supplied free electrons */ 01233 gv.TotalEden ); 01234 01235 /* average charge per grain in electrons */ 01236 for( nd=0; nd < gv.nBin; ++nd ) 01237 { 01238 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ ); 01239 } 01240 fprintf( punch.ipPnunit[ipPun], "\n" ); 01241 } 01242 } 01243 01244 else if( strcmp(punch.chPunch[ipPun],"DUSH") == 0 ) 01245 { 01246 /* grain heating */ 01247 if( strcmp(chTime,"LAST") != 0 ) 01248 { 01249 /* used to print header exactly one time */ 01250 static bool lgMustPrtHeader = true; 01251 /* punch grain charge, but do labels first if this is first zone */ 01252 if( lgMustPrtHeader ) 01253 { 01254 /* first print string giving grain id */ 01255 fprintf( punch.ipPnunit[ipPun], "#Depth" ); 01256 for( nd=0; nd < gv.nBin; ++nd ) 01257 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); 01258 fprintf( punch.ipPnunit[ipPun], "\n" ); 01259 01260 /* now print grain radius converting from cm to microns */ 01261 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" ); 01262 for( nd=0; nd < gv.nBin; ++nd ) 01263 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); 01264 fprintf( punch.ipPnunit[ipPun], "\n" ); 01265 01266 /* don't need to do this, ever again */ 01267 lgMustPrtHeader = false; 01268 } 01269 fprintf( punch.ipPnunit[ipPun], " %.5e", 01270 radius.depth_mid_zone ); 01271 /* grain heating */ 01272 for( nd=0; nd < gv.nBin; ++nd ) 01273 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl ); 01274 fprintf( punch.ipPnunit[ipPun], "\n" ); 01275 } 01276 } 01277 01278 else if( strcmp(punch.chPunch[ipPun],"DUSV") == 0 ) 01279 { 01280 /* grain drift velocities */ 01281 if( strcmp(chTime,"LAST") != 0 ) 01282 { 01283 /* used to print header exactly one time */ 01284 static bool lgMustPrtHeader = true; 01285 /* punch grain velocity, but do labels first if this is first zone */ 01286 if( lgMustPrtHeader ) 01287 { 01288 /* first print string giving grain id */ 01289 fprintf( punch.ipPnunit[ipPun], "#Depth" ); 01290 for( nd=0; nd < gv.nBin; ++nd ) 01291 fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab ); 01292 fprintf( punch.ipPnunit[ipPun], "\n" ); 01293 01294 /* now print grain radius converting from cm to microns */ 01295 fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" ); 01296 for( nd=0; nd < gv.nBin; ++nd ) 01297 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 ); 01298 fprintf( punch.ipPnunit[ipPun], "\n" ); 01299 01300 /* don't need to do this, ever again */ 01301 lgMustPrtHeader = false; 01302 } 01303 fprintf( punch.ipPnunit[ipPun], " %.5e", 01304 radius.depth_mid_zone ); 01305 /* grain drift velocity in km/s */ 01306 for( nd=0; nd < gv.nBin; ++nd ) 01307 fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 ); 01308 fprintf( punch.ipPnunit[ipPun], "\n" ); 01309 } 01310 } 01311 01312 /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */ 01313 else if( strcmp(punch.chPunch[ipPun],"DUSQ") == 0 ) 01314 { 01315 /* punch grain Qs */ 01316 if( strcmp(chTime,"LAST") == 0 ) 01317 { 01318 for( j=0; j < rfield.nflux; j++ ) 01319 { 01320 fprintf( punch.ipPnunit[ipPun], "%11.4e", 01321 rfield.anu[j] ); 01322 for( nd=0; nd < gv.nBin; nd++ ) 01323 { 01324 fprintf( punch.ipPnunit[ipPun], "%9.1e%9.1e", 01325 gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea, 01326 gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea ); 01327 } 01328 fprintf( punch.ipPnunit[ipPun], "\n" ); 01329 } 01330 } 01331 } 01332 01333 else if( strcmp(punch.chPunch[ipPun],"ELEM") == 0 ) 01334 { 01335 if( strcmp(chTime,"LAST") != 0 ) 01336 { 01337 realnum renorm = 1.f; 01338 01339 /* this is the index for the atomic number on the physical scale */ 01340 /* >>chng 04 nov 23, use c scale throughout */ 01341 nelem = (long int)punch.punarg[ipPun][0]; 01342 ASSERT( nelem >= ipHYDROGEN ); 01343 01344 /* don't do this if element is not turned on */ 01345 if( dense.lgElmtOn[nelem] ) 01346 { 01347 /* >>chng 04 nov 23, add density option, leave as cm-3 01348 * default is still norm to total of that element */ 01349 if( punch.punarg[ipPun][1] == 0 ) 01350 renorm = dense.gas_phase[nelem]; 01351 01352 fprintf( punch.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone ); 01353 01354 for( j=0; j <= (nelem + 1); ++j) 01355 { 01356 fprintf( punch.ipPnunit[ipPun], "\t%.2e", 01357 dense.xIonDense[nelem][j]/renorm ); 01358 } 01359 if( nelem==ipHYDROGEN ) 01360 { 01361 /* H2 */ 01362 fprintf( punch.ipPnunit[ipPun], "\t%.2e", 01363 hmi.H2_total/renorm ); 01364 } 01365 /* >>chng 04 nov 23 add C and O fine structure pops */ 01366 else if( nelem==ipCARBON ) 01367 { 01368 fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e", 01369 colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm); 01370 fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e", 01371 colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm); 01372 fprintf( punch.ipPnunit[ipPun], "\t%.2e", 01373 findspecies("CO")->hevmol/renorm ); 01374 } 01375 else if( nelem==ipOXYGEN ) 01376 { 01377 fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e", 01378 colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm); 01379 } 01380 fprintf( punch.ipPnunit[ipPun], "\n" ); 01381 } 01382 } 01383 } 01384 01385 else if( strcmp(punch.chPunch[ipPun],"RECA") == 0 ) 01386 { 01387 /* this will create table for AGN3 then exit, 01388 * routine is in makerecom.c */ 01389 ion_recombAGN( punch.ipPnunit[ipPun] ); 01390 cdEXIT(EXIT_FAILURE); 01391 } 01392 01393 else if( strcmp(punch.chPunch[ipPun],"RECE") == 0 ) 01394 { 01395 /* punch recombination efficiencies, 01396 * option turned on with the "punch recombination efficiencies" command 01397 * output for the punch recombination coefficients command is actually 01398 * produced by a series of routines, as they generate the recombination 01399 * coefficients. these include 01400 * dielsupres, helium, hydrorecom, iibod, and makerecomb*/ 01401 fprintf( punch.ipPnunit[ipPun], 01402 "%12.4e %12.4e %12.4e %12.4e\n", 01403 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][0][ipRecRad], 01404 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][0][ipRecNetEsc] , 01405 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecRad], 01406 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecNetEsc]); 01407 } 01408 01409 else if( strcmp(punch.chPunch[ipPun],"FEco") == 0 ) 01410 { 01411 /* FeII continuum, check that FeII is turned on */ 01412 if( strcmp(chTime,"LAST") == 0 ) 01413 { 01414 for( j=0; j < nFeIIConBins; j++ ) 01415 { 01416 /* emission from large FeII atom, integrated over continuum intervals 01417 * units are always in units, erg cm-2 s-1 as in intensity case, 01418 * even if luminosity case is considered 01419 * [j][0] is flux, [j][1] and [j][2] are upper and 01420 * lower bounds in Angstroms. 01421 * these are set in FeIIZero */ 01424 fprintf( punch.ipPnunit[ipPun], "%.2f\t%e \n", 01425 (FeII_Cont[j][1]+FeII_Cont[j][2])/2., 01426 FeII_Cont[j][0] ); 01427 } 01428 } 01429 } 01430 01431 /* punch column densities */ 01432 else if( strcmp(punch.chPunch[ipPun],"FEcl") == 0 ) 01433 { 01434 if( strcmp(chTime,"LAST") == 0 ) 01435 { 01436 /* punch FeII level energies and stat weights, followed by column density */ 01437 FeIIPunchColden( punch.ipPnunit[ipPun] ); 01438 } 01439 } 01440 01441 else if( strcmp(punch.chPunch[ipPun],"FE2l") == 0 ) 01442 { 01443 if( strcmp(chTime,"LAST") == 0 ) 01444 { 01445 /* punch FeII level energies and stat weights */ 01446 FeIIPunchLevels( punch.ipPnunit[ipPun] ); 01447 } 01448 } 01449 01450 else if( strcmp(punch.chPunch[ipPun],"FEli") == 0 ) 01451 { 01452 if( strcmp(chTime,"LAST") == 0 ) 01453 { 01454 /* punch line intensities, routine is in atom_feii.c */ 01455 FeIIPunchLines( punch.ipPnunit[ipPun] ); 01456 } 01457 } 01458 01459 else if( strcmp(punch.chPunch[ipPun],"FE2o") == 0 ) 01460 { 01461 if( strcmp(chTime,"LAST") == 0 ) 01462 { 01463 /* punch line optical depths, routine is in atom_feii.c */ 01464 FeIIPunchOpticalDepth( punch.ipPnunit[ipPun] ); 01465 } 01466 } 01467 01468 else if( strcmp(punch.chPunch[ipPun],"FRED") == 0 ) 01469 { 01470 /* set with punch Fred command, this punches some stuff from 01471 * Fred Hamann's dynamics project */ 01472 if( strcmp(chTime,"LAST") != 0 ) 01473 { 01474 /* Fred's list */ 01475 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\t%.3e\t%.3e\t%.3e" 01476 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e" 01477 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e" 01478 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e" 01479 "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 01480 radius.Radius, radius.depth ,wind.windv/1e5, 01481 dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te, 01482 wind.AccelLine , wind.AccelCont , 01483 wind.fmul , 01484 mean.xIonMeans[0][ipHYDROGEN][0] , mean.xIonMeans[0][ipHYDROGEN][1] , 01485 mean.xIonMeans[0][ipHELIUM][0] , mean.xIonMeans[0][ipHELIUM][1] , 01486 mean.xIonMeans[0][ipHELIUM][2] , 01487 mean.xIonMeans[0][ipCARBON][1] , mean.xIonMeans[0][ipCARBON][2] , 01488 mean.xIonMeans[0][ipCARBON][3] , 01489 mean.xIonMeans[0][ipOXYGEN][0] , mean.xIonMeans[0][ipOXYGEN][1] , 01490 mean.xIonMeans[0][ipOXYGEN][2] , mean.xIonMeans[0][ipOXYGEN][3] , 01491 mean.xIonMeans[0][ipOXYGEN][4] , mean.xIonMeans[0][ipOXYGEN][5] , 01492 mean.xIonMeans[0][ipOXYGEN][6] , mean.xIonMeans[0][ipOXYGEN][7] , 01493 dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] , 01494 dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] , 01495 dense.xIonDense[ipHELIUM][2] , 01496 dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] , 01497 dense.xIonDense[ipCARBON][3] , 01498 dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] , 01499 dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] , 01500 dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] , 01501 dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] , 01502 mean.xIonMeans[0][ipMAGNESIUM][1] , dense.xIonDense[ipMAGNESIUM][1]); 01503 } 01504 } 01505 01506 else if( strcmp(punch.chPunch[ipPun],"FE2d") == 0 ) 01507 { 01508 /* punch some departure coefficients for large FeII atom */ 01509 if( strcmp(chTime,"LAST") != 0 ) 01510 FeIIPunDepart(punch.ipPnunit[ipPun],false); 01511 } 01512 01513 else if( strcmp(punch.chPunch[ipPun],"FE2D") == 0 ) 01514 { 01515 /* punch all departure coefficients for large FeII atom */ 01516 if( strcmp(chTime,"LAST") != 0 ) 01517 FeIIPunDepart(punch.ipPnunit[ipPun],true); 01518 } 01519 01520 else if( strcmp(punch.chPunch[ipPun],"FE2p") == 0 ) 01521 { 01522 bool lgFlag = false; 01523 if( punch.punarg[ipPun][2] ) 01524 lgFlag = true; 01525 /* punch small subset of level populations for large FeII atom */ 01526 if( strcmp(chTime,"LAST") != 0 ) 01527 FeIIPunPop(punch.ipPnunit[ipPun],false,0,0, 01528 lgFlag); 01529 } 01530 01531 else if( strcmp(punch.chPunch[ipPun],"FE2P") == 0 ) 01532 { 01533 bool lgFlag = false; 01534 if( punch.punarg[ipPun][2] ) 01535 lgFlag = true; 01536 /* punch range of level populations for large FeII atom */ 01537 if( strcmp(chTime,"LAST") != 0 ) 01538 FeIIPunPop(punch.ipPnunit[ipPun], 01539 true, 01540 (long int)punch.punarg[ipPun][0], 01541 (long int)punch.punarg[ipPun][1], 01542 lgFlag); 01543 } 01544 01545 /* punch spectra in fits format */ 01546 else if( strcmp(punch.chPunch[ipPun],"FITS") == 0 ) 01547 { 01548 if( strcmp(chTime,"LAST") == 0 ) 01549 { 01550 punchFITSfile( punch.ipPnunit[ipPun], NUM_OUTPUT_TYPES ); 01551 } 01552 } 01553 /* punch gammas (but without element) */ 01554 else if( strcmp(punch.chPunch[ipPun],"GAMt") == 0 ) 01555 { 01556 if( strcmp(chTime,"LAST") != 0 ) 01557 { 01558 long ns; 01559 /* punch photoionization rates, with the PUNCH GAMMAS command */ 01560 for( nelem=0; nelem < LIMELM; nelem++ ) 01561 { 01562 for( ion=0; ion <= nelem; ion++ ) 01563 { 01564 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 01565 { 01566 fprintf( punch.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e", 01567 nelem+1, ion+1, ns+1, 01568 ionbal.PhotoRate_Shell[nelem][ion][ns][0], 01569 ionbal.PhotoRate_Shell[nelem][ion][ns][1] , 01570 ionbal.PhotoRate_Shell[nelem][ion][ns][2] ); 01571 01572 for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ ) 01573 { 01574 fprintf( punch.ipPnunit[ipPun], "%5.2f", 01575 t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) ); 01576 } 01577 fprintf( punch.ipPnunit[ipPun], "\n" ); 01578 } 01579 } 01580 } 01581 } 01582 } 01583 01584 /* punch gammas element, ion */ 01585 else if( strcmp(punch.chPunch[ipPun],"GAMe") == 0 ) 01586 { 01587 if( strcmp(chTime,"LAST") != 0 ) 01588 { 01589 int ns; 01590 nelem = (long)punch.punarg[ipPun][0]; 01591 ion = (long)punch.punarg[ipPun][1]; 01592 /* valence shell */ 01593 ns = Heavy.nsShells[nelem][ion]-1; 01594 /* show what some of the ionization sources are */ 01595 GammaPrt( 01596 opac.ipElement[nelem][ion][ns][0] , 01597 opac.ipElement[nelem][ion][ns][1] , 01598 opac.ipElement[nelem][ion][ns][2] , 01599 punch.ipPnunit[ipPun], 01600 ionbal.PhotoRate_Shell[nelem][ion][ns][0] , 01601 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 ); 01602 } 01603 } 01604 01605 else if( strcmp(punch.chPunch[ipPun],"GAUN") == 0 ) 01606 { 01607 /* punch gaunt factors */ 01608 if( strcmp(chTime,"LAST") != 0 ) 01609 PunchGaunts(punch.ipPnunit[ipPun]); 01610 } 01611 01612 /* punch parameters of grid calculation */ 01613 else if( strcmp(punch.chPunch[ipPun],"GRID") == 0 ) 01614 { 01615 /* one time punch of all grid parameters after grid is complete */ 01616 if( (strcmp(chTime,"LAST") == 0 ) && grid.lgGridDone ) 01617 { 01618 /* start of line gives abort and warning summary */ 01619 fprintf(punch.ipPnunit[ipPun], "#Abort?\tWarnings?" ); 01620 /* print start of each variable command line */ 01621 for( i=0; i< grid.nintparm; i++ ) 01622 { 01623 char chStr[10]; 01624 strncpy( chStr , optimize.chVarFmt[i] , 9 ); 01625 /* make sure this small bit of string is terminated */ 01626 chStr[9] = '\0'; 01627 fprintf(punch.ipPnunit[ipPun], "\t%s", chStr ); 01628 } 01629 fprintf(punch.ipPnunit[ipPun], "\n" ); 01630 for( i=0; i<grid.totNumModels; i++ ) 01631 { 01632 /* abort / warning summary for this sim */ 01633 fprintf(punch.ipPnunit[ipPun], "%c\t%c", 01634 TorF(grid.lgAbort[i]) , 01635 TorF(grid.lgWarn[i]) ); 01636 /* the grid parameters */ 01637 for( j=0; j< grid.nintparm; j++ ) 01638 { 01639 fprintf(punch.ipPnunit[ipPun], "\t%f", grid.interpParameters[i][j] ); 01640 } 01641 fprintf(punch.ipPnunit[ipPun], "\n" ); 01642 } 01643 } 01644 } 01645 01646 else if( strcmp(punch.chPunch[ipPun],"HISp") == 0 ) 01647 { 01648 /* punch pressure history of current zone */ 01649 if( strcmp(chTime,"LAST") != 0 ) 01650 { 01651 /* note if pressure convergence failure occurred in history that follows */ 01652 if( !conv.lgConvPres ) 01653 { 01654 fprintf( punch.ipPnunit[ipPun], 01655 "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n", 01656 iteration , nzone ); 01657 } 01658 /* note if temperature convergence failure occurred in history that follows */ 01659 if( !conv.lgConvTemp ) 01660 { 01661 fprintf( punch.ipPnunit[ipPun], 01662 "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n", 01663 iteration , nzone ); 01664 } 01665 for(i=0; i<conv.hist_pres_npres; ++i ) 01666 { 01667 /* punch history of density - pressure, with correct pressure */ 01668 fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n", 01669 iteration, 01670 nzone, 01671 conv.hist_pres_density[i], 01672 conv.hist_pres_current[i], 01673 conv.hist_pres_correct[i]); 01674 } 01675 } 01676 } 01677 01678 else if( strcmp(punch.chPunch[ipPun],"HISt") == 0 ) 01679 { 01680 /* punch temperature history of current zone */ 01681 if( strcmp(chTime,"LAST") != 0 ) 01682 { 01683 /* note if pressure convergence failure occurred in history that follows */ 01684 if( !conv.lgConvPres ) 01685 { 01686 fprintf( punch.ipPnunit[ipPun], 01687 "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n", 01688 iteration , nzone ); 01689 } 01690 /* note if temperature convergence failure occurred in history that follows */ 01691 if( !conv.lgConvTemp ) 01692 { 01693 fprintf( punch.ipPnunit[ipPun], 01694 "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n", 01695 iteration , nzone ); 01696 } 01697 for(i=0; i<conv.hist_temp_ntemp; ++i ) 01698 { 01699 /* punch history of density - pressure, with correct pressure */ 01700 fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n", 01701 iteration, 01702 nzone, 01703 conv.hist_temp_temp[i], 01704 conv.hist_temp_heat[i], 01705 conv.hist_temp_cool[i]); 01706 } 01707 } 01708 } 01709 01710 else if( strncmp(punch.chPunch[ipPun],"H2",2) == 0 ) 01711 { 01712 /* all punch info on large H2 molecule include H2 PDR pdr */ 01713 H2_PunchDo( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] , chTime, ipPun ); 01714 } 01715 01716 else if( strcmp(punch.chPunch[ipPun],"HEAT") == 0 ) 01717 { 01718 /* punch heating, routine in file of same name */ 01719 if( strcmp(chTime,"LAST") != 0 ) 01720 HeatPunch(punch.ipPnunit[ipPun]); 01721 } 01722 01723 else if( strncmp(punch.chPunch[ipPun],"HE",2) == 0 ) 01724 { 01725 /* various punch helium commands */ 01726 /* punch helium line wavelengths */ 01727 if( strcmp(punch.chPunch[ipPun] , "HELW") == 0 ) 01728 { 01729 if( strcmp(chTime,"LAST") == 0 ) 01730 { 01731 /* punch helium & he-like wavelengths, first header */ 01732 fprintf( punch.ipPnunit[ipPun], 01733 "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S" 01734 "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" ); 01735 fprintf( punch.ipPnunit[ipPun], "\n" ); 01736 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 01737 { 01738 /* print element name, nuclear charge */ 01739 fprintf( punch.ipPnunit[ipPun], "%li\t%s", 01740 nelem+1 , elementnames.chElementSym[nelem] ); 01741 /*prt_wl print floating wavelength in Angstroms, in output format */ 01742 fprintf( punch.ipPnunit[ipPun], "\t" ); 01743 prt_wl( punch.ipPnunit[ipPun] , 01744 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].WLAng ); 01745 fprintf( punch.ipPnunit[ipPun], "\t" ); 01746 prt_wl( punch.ipPnunit[ipPun] , 01747 Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe1s1S].WLAng ); 01748 fprintf( punch.ipPnunit[ipPun], "\t" ); 01749 prt_wl( punch.ipPnunit[ipPun] , 01750 Transitions[ipHE_LIKE][nelem][ipHe2p3P2][ipHe1s1S].WLAng ); 01751 fprintf( punch.ipPnunit[ipPun], "\t" ); 01752 prt_wl( punch.ipPnunit[ipPun] , 01753 Transitions[ipHE_LIKE][nelem][ipHe2s3S][ipHe1s1S].WLAng ); 01754 fprintf( punch.ipPnunit[ipPun], "\t" ); 01755 prt_wl( punch.ipPnunit[ipPun] , 01756 Transitions[ipHE_LIKE][nelem][ipHe2p3P2][ipHe2s3S].WLAng ); 01757 fprintf( punch.ipPnunit[ipPun], "\t" ); 01758 prt_wl( punch.ipPnunit[ipPun] , 01759 Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe2s3S].WLAng ); 01760 fprintf( punch.ipPnunit[ipPun], "\t" ); 01761 prt_wl( punch.ipPnunit[ipPun] , 01762 Transitions[ipHE_LIKE][nelem][ipHe2p3P0][ipHe2s3S].WLAng ); 01763 fprintf( punch.ipPnunit[ipPun], "\n"); 01764 } 01765 } 01766 } 01767 else 01768 TotalInsanity(); 01769 } 01770 01771 /* punch hummer, results needed for Lya transport, to feed into David's routine */ 01772 else if( strcmp(punch.chPunch[ipPun],"HUMM") == 0 ) 01773 { 01774 eps = Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul/ 01775 (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL *dense.eden); 01776 fprintf( punch.ipPnunit[ipPun], 01777 " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n", 01778 radius.depth_mid_zone, 01779 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn, 01780 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1], 01781 StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*dense.xIonDense[ipHYDROGEN][1], 01782 phycon.te, Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->damp, eps ); 01783 } 01784 01785 else if( strncmp( punch.chPunch[ipPun] , "HYD", 3 ) == 0 ) 01786 { 01787 /* various punch hydrogen commands */ 01788 if( strcmp(punch.chPunch[ipPun],"HYDc") == 0 ) 01789 { 01790 if( strcmp(chTime,"LAST") != 0 ) 01791 { 01792 /* punch hydrogen physical conditions */ 01793 fprintf( punch.ipPnunit[ipPun], 01794 " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 01795 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden, 01796 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN], 01797 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN], 01798 hmi.H2_total/dense.gas_phase[ipHYDROGEN], 01799 hmi.Hmolec[ipMH2p]/dense.gas_phase[ipHYDROGEN], 01800 hmi.Hmolec[ipMH3p]/dense.gas_phase[ipHYDROGEN], 01801 hmi.Hmolec[ipMHm]/dense.gas_phase[ipHYDROGEN] ); 01802 } 01803 } 01804 01805 else if( strcmp(punch.chPunch[ipPun],"HYDi") == 0 ) 01806 { 01807 if( strcmp(chTime,"LAST") != 0 ) 01808 { 01809 /* punch hydrogen ionization 01810 * this will be total decays to ground */ 01811 RateInter = 0.; 01812 stage = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][0]*dense.eden*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop; 01813 fref = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop; 01814 fout = StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop; 01815 /* 06 aug 28, from numLevels_max to _local. */ 01816 for( ion=ipH2s; ion < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ion++ ) 01817 { 01818 /* this is total decays to ground */ 01819 RateInter += 01820 Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Aul* 01821 (Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pesc + 01822 Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pelec_esc + 01823 Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pdest); 01824 /* total photo from all levels */ 01825 fref += iso.gamnc[ipH_LIKE][ipHYDROGEN][ion]*StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop; 01826 /* total col ion from all levels */ 01827 stage += iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ion]*dense.eden* 01828 StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop; 01829 fout += StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop; 01830 } 01831 fprintf( punch.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e", 01832 nzone, 01833 /* photo and collision ion rates have units s-1 */ 01834 iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s], 01835 iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s]* dense.EdenHCorr, 01836 ionbal.RateRecomTot[ipHYDROGEN][0] ); 01837 01838 fprintf( punch.ipPnunit[ipPun], "\t%.2e", 01839 iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN] ); 01840 01841 fprintf( punch.ipPnunit[ipPun], 01842 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 01843 dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0], 01844 iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/(ionbal.RateRecomTot[ipHYDROGEN][0]), 01845 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][1][ipRecEsc], 01846 RateInter, 01847 fref/MAX2(1e-37,fout), 01848 stage/MAX2(1e-37,fout), 01849 /* simple H+ */ 01850 safe_div( iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][0], dense.eden*dense.xIonDense[ipHYDROGEN][1] ), 01851 secondaries.csupra[ipHYDROGEN][0]); 01852 01853 GammaPrt(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0],rfield.nflux,iso.ipOpac[ipH_LIKE][ipHYDROGEN][ipH1s], 01854 punch.ipPnunit[ipPun],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]* 01855 0.05); 01856 } 01857 } 01858 01859 else if( strcmp(punch.chPunch[ipPun],"HYDp") == 0 ) 01860 { 01861 if( strcmp(chTime,"LAST") != 0 ) 01862 { 01863 /* punch hydrogen populations 01864 * first give total atom and ion density [cm-3]*/ 01865 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e", 01866 radius.depth_mid_zone, 01867 dense.xIonDense[ipHYDROGEN][0], 01868 dense.xIonDense[ipHYDROGEN][1] ); 01869 01870 /* next give state-specific densities [cm-3] */ 01871 for( j=ipH1s; j < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]-1; j++ ) 01872 { 01873 fprintf( punch.ipPnunit[ipPun], "\t%.2e", 01874 StatesElem[ipH_LIKE][ipHYDROGEN][j].Pop*dense.xIonDense[ipHYDROGEN][1] ); 01875 } 01876 fprintf( punch.ipPnunit[ipPun], "\n" ); 01877 } 01878 } 01879 01880 else if( strcmp(punch.chPunch[ipPun],"HYDl") == 0 ) 01881 { 01882 if( strcmp(chTime,"LAST") == 0 ) 01883 { 01884 /* punch hydrogen line 01885 * gives intensities and optical depths */ 01888 for( ipHi=4; ipHi<iso.numLevels_local[ipH_LIKE][ipHYDROGEN]-1; ++ipHi ) 01889 { 01890 for( ipLo=ipHi-5; ipLo<ipHi; ++ipLo ) 01891 { 01892 if( ipLo<0 ) 01893 continue; 01894 fprintf(punch.ipPnunit[ipPun], "%li\t%li\t%.4e\t%.2e\n", 01895 ipHi, ipLo, 01896 /* wl in cm */ 01897 1./Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].EnergyWN, 01898 Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Emis->TauIn ); 01899 } 01900 } 01901 } 01902 } 01903 01904 /* punch hydrogen Lya - some details about Lya */ 01905 else if( strcmp(punch.chPunch[ipPun],"HYDL") == 0 ) 01906 { 01907 if( strcmp(chTime,"LAST") != 0 ) 01908 { 01909 /* the population ratio for Lya */ 01910 double popul = StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop/SDIV(StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop); 01911 /* the excitation temperature of Lya */ 01912 texc = TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ); 01913 fprintf( punch.ipPnunit[ipPun], 01914 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 01915 radius.depth_mid_zone, 01916 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn, 01917 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauTot, 01918 popul, 01919 texc, 01920 phycon.te, 01921 texc/phycon.te , 01922 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pesc, 01923 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest, 01924 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->pump, 01925 opac.opacity_abs[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1], 01926 opac.albedo[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ); 01927 } 01928 } 01929 01930 else if( strcmp(punch.chPunch[ipPun],"HYDr") == 0 ) 01931 { 01932 /* punch hydrogen recc - recombination cooling for AGN3 */ 01933 TempChange(2500.f, false); 01934 while( phycon.te <= 20000. ) 01935 { 01936 double r1; 01937 double ThinCoolingCaseB; 01938 01939 r1 = HydroRecCool(1,0); 01940 ThinCoolingCaseB = pow(10.,((-25.859117 + 01941 0.16229407*phycon.telogn[0] + 01942 0.34912863*phycon.telogn[1] - 01943 0.10615964*phycon.telogn[2])/(1. + 01944 0.050866793*phycon.telogn[0] - 01945 0.014118924*phycon.telogn[1] + 01946 0.0044980897*phycon.telogn[2] + 01947 6.0969594e-5*phycon.telogn[3])))/phycon.te; 01948 01949 fprintf( punch.ipPnunit[ipPun], " %10.2e\t", 01950 phycon.te); 01951 fprintf( punch.ipPnunit[ipPun], " %10.2e\t", 01952 (r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) ); 01953 01954 fprintf( punch.ipPnunit[ipPun], " %10.2e\t", 01955 r1/(BOLTZMANN*phycon.te)); 01956 01957 fprintf( punch.ipPnunit[ipPun], " %10.2e\n", 01958 ThinCoolingCaseB/(BOLTZMANN*phycon.te)); 01959 01960 TempChange(phycon.te *2.f , false); 01961 } 01962 /* must exit since we have disturbed the solution */ 01963 fprintf(ioQQQ , "punch agn now exits since solution is disturbed.\n"); 01964 cdEXIT( EXIT_SUCCESS ); 01965 } 01966 else 01967 TotalInsanity(); 01968 } 01969 01970 else if( strcmp(punch.chPunch[ipPun],"IONI") == 0 ) 01971 { 01972 if( strcmp(chTime,"LAST") == 0 ) 01973 { 01974 /* punch mean ionization distribution */ 01975 PrtMeanIon( 'i', false , punch.ipPnunit[ipPun] ); 01976 } 01977 } 01978 01979 /* punch ionization rates */ 01980 else if( strcmp(punch.chPunch[ipPun],"IONR") == 0 ) 01981 { 01982 if( strcmp(chTime,"LAST") != 0 ) 01983 { 01984 /* this is element number */ 01985 nelem = (long)punch.punarg[ipPun][0]; 01986 fprintf( punch.ipPnunit[ipPun], 01987 "%.5e\t%.4e\t%.4e", 01988 radius.depth_mid_zone, 01989 dense.eden , 01990 dynamics.Rate); 01991 /* >>chng 04 oct 15, from nelem+2 to nelem+1 - array over read - 01992 * caught by PnH */ 01993 for( ion=0; ion<nelem+1; ++ion ) 01994 { 01995 fprintf( punch.ipPnunit[ipPun], 01996 "\t%.4e\t%.4e\t%.4e\t%.4e", 01997 dense.xIonDense[nelem][ion] , 01998 ionbal.RateIonizTot[nelem][ion] , 01999 ionbal.RateRecomTot[nelem][ion] , 02000 dynamics.Source[nelem][ion] ); 02001 } 02002 fprintf( punch.ipPnunit[ipPun], "\n"); 02003 } 02004 } 02005 02006 else if( strcmp(punch.chPunch[ipPun]," IP ") == 0 ) 02007 { 02008 if( strcmp(chTime,"LAST") == 0 ) 02009 { 02010 /* punch valence shell ip's */ 02011 for( nelem=0; nelem < LIMELM; nelem++ ) 02012 { 02013 int ion_big; 02014 double energy; 02015 02016 /* this is the largest number of ion stages per line */ 02017 const int NELEM_LINE = 10; 02018 /* this loop in case all ions do not fit across page */ 02019 for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE ) 02020 { 02021 int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem); 02022 02023 /* new line then element name */ 02024 fprintf( punch.ipPnunit[ipPun], 02025 "\n%2.2s", elementnames.chElementSym[nelem]); 02026 02027 /* print ion stages across line */ 02028 for( ion=ion_big; ion <= ion_limit; ++ion ) 02029 { 02030 fprintf( punch.ipPnunit[ipPun], "\t%4ld", ion+1 ); 02031 } 02032 fprintf( punch.ipPnunit[ipPun], "\n" ); 02033 02034 /* this loop is over all shells */ 02035 ASSERT( ion_limit < LIMELM ); 02036 /* upper limit is number of shells in atom */ 02037 for( ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ ) 02038 { 02039 02040 /* print shell label */ 02041 fprintf( punch.ipPnunit[ipPun], "%2.2s", Heavy.chShell[ips]); 02042 02043 /* loop over possible ions */ 02044 for( ion=ion_big; ion<=ion_limit; ++ion ) 02045 { 02046 02047 /* does this subshell exist for this ion? break if it does not*/ 02048 /*if( Heavy.nsShells[nelem][ion]<Heavy.nsShells[nelem][0] )*/ 02049 if( ips >= Heavy.nsShells[nelem][ion] ) 02050 break; 02051 02052 /* array elements are shell, numb of electrons, element, 0 */ 02053 energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0); 02054 02055 /* now print threshold with correct format */ 02056 if( energy < 10. ) 02057 { 02058 fprintf( punch.ipPnunit[ipPun], "\t%6.3f", energy ); 02059 } 02060 else if( energy < 100. ) 02061 { 02062 fprintf( punch.ipPnunit[ipPun], "\t%6.2f", energy ); 02063 } 02064 else if( energy < 1000. ) 02065 { 02066 fprintf( punch.ipPnunit[ipPun], "\t%6.1f", energy ); 02067 } 02068 else 02069 { 02070 fprintf( punch.ipPnunit[ipPun], "\t%6ld", (long)(energy) ); 02071 } 02072 } 02073 02074 /* put cs at end of long line */ 02075 fprintf( punch.ipPnunit[ipPun], "\n" ); 02076 } 02077 } 02078 } 02079 } 02080 } 02081 02082 else if( strcmp(punch.chPunch[ipPun],"LINC") == 0 ) 02083 { 02084 /* punch line cumulative */ 02085 /* lgOK not used, placdholder */ 02086 if( strcmp(chTime,"LAST") != 0 ) 02087 { 02088 /* not used for anything here, but should not pass baloney*/ 02089 lgOK = true; 02090 punch_line(punch.ipPnunit[ipPun],"PUNC",lgOK,chDummy); 02091 } 02092 } 02093 02094 else if( strcmp(punch.chPunch[ipPun],"LIND") == 0 ) 02095 { 02096 /* punch line data, then stop */ 02097 PunchLineData(punch.ipPnunit[ipPun]); 02098 } 02099 02100 else if( strcmp(punch.chPunch[ipPun],"LINL") == 0 ) 02101 { 02102 /* punch line labels */ 02103 bool lgPrintAll=false; 02104 /* LONG keyword on punch line labels command sets this to 1 */ 02105 if( punch.punarg[ipPun][0]>0. ) 02106 lgPrintAll = true; 02107 prt_LineLabels(punch.ipPnunit[ipPun] , lgPrintAll ); 02108 } 02109 02110 else if( strcmp(punch.chPunch[ipPun],"LINO") == 0 ) 02111 { 02112 if( strcmp(chTime,"LAST") == 0 ) 02113 { 02114 /* punch line optical depths, routine is below, file static */ 02115 PunchLineStuff(punch.ipPnunit[ipPun],"optical" , punch.punarg[ipPun][0]); 02116 } 02117 } 02118 02119 else if( strcmp(punch.chPunch[ipPun],"LINP") == 0 ) 02120 { 02121 if( strcmp(chTime,"LAST") != 0 ) 02122 { 02123 static bool lgFirst=true; 02124 /* punch line populations, need to do this twice if very first 02125 * call since first call to PunchLineStuff generates atomic parameters 02126 * rather than level pops, routine is below, file static */ 02127 PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[ipPun][0]); 02128 if( lgFirst ) 02129 { 02130 lgFirst = false; 02131 PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[ipPun][0]); 02132 } 02133 } 02134 } 02135 02136 else if( strcmp(punch.chPunch[ipPun],"LINS") == 0 ) 02137 { 02138 /* punch line emissivity */ 02139 if( strcmp(chTime,"LAST") != 0 ) 02140 { 02141 /* not actually used, but should not pass baloney*/ 02142 lgOK = true; 02143 punch_line(punch.ipPnunit[ipPun],"PUNS",lgOK,chDummy); 02144 } 02145 } 02146 02147 else if( strcmp(punch.chPunch[ipPun],"LINR") == 0 ) 02148 { 02149 /* punch line RT */ 02150 if( strcmp(chTime,"LAST") != 0 ) 02151 { 02152 Punch_Line_RT( punch.ipPnunit[ipPun] , "PUNC" ); 02153 } 02154 } 02155 02156 else if( strcmp(punch.chPunch[ipPun],"LINA") == 0 ) 02157 { 02158 /* punch line array */ 02159 if( strcmp(chTime,"LAST") == 0 ) 02160 { 02161 /* punch out all lines with energies */ 02162 for( j=0; j < LineSave.nsum; j++ ) 02163 { 02164 if( LineSv[j].wavelength > 0. && 02165 LineSv[j].sumlin[LineSave.lgLineEmergent] > 0. ) 02166 { 02167 /* line energy, in units set with units option */ 02168 fprintf( punch.ipPnunit[ipPun], "%12.5e", 02169 AnuUnit((realnum)RYDLAM/LineSv[j].wavelength) ); 02170 /* line label */ 02171 fprintf( punch.ipPnunit[ipPun], "\t%4.4s ", 02172 LineSv[j].chALab ); 02173 /* wavelength */ 02174 prt_wl( punch.ipPnunit[ipPun], LineSv[j].wavelength ); 02175 /* intrinsic intensity */ 02176 fprintf( punch.ipPnunit[ipPun], "\t%8.3f", 02177 log10(SDIV(LineSv[j].sumlin[0]) ) + radius.Conv2PrtInten ); 02178 /* emergent line intensity, r recombination */ 02179 fprintf( punch.ipPnunit[ipPun], "\t%8.3f", 02180 log10(SDIV(LineSv[j].sumlin[1]) ) + radius.Conv2PrtInten ); 02181 /* type of line, i for info, etc */ 02182 fprintf( punch.ipPnunit[ipPun], " \t%c\n", 02183 LineSv[j].chSumTyp); 02184 } 02185 } 02186 } 02187 } 02188 02189 else if( strcmp(punch.chPunch[ipPun],"LINI") == 0 ) 02190 { 02191 if( strcmp(chTime,"LAST") == 0 && 02192 (nzone/punch.LinEvery)*punch.LinEvery != nzone ) 02193 { 02194 /* this is the last zone 02195 * punch line intensities - but do not do last zone twice */ 02196 PunLineIntensity(punch.ipPnunit[ipPun]); 02197 } 02198 else if( strcmp(chTime,"LAST") != 0 ) 02199 { 02200 /* following so we only punch first zone if LinEvery reset */ 02201 if( (punch.lgLinEvery && nzone == 1) || 02202 (nzone/punch.LinEvery)*punch.LinEvery == 02203 nzone ) 02204 { 02205 /* this is middle of calculation 02206 * punch line intensities */ 02207 PunLineIntensity(punch.ipPnunit[ipPun]); 02208 } 02209 } 02210 } 02211 02212 else if( strcmp( punch.chPunch[ipPun],"LEIL") == 0) 02213 { 02214 /* some line intensities for the Leiden PDR, 02215 * but only do this when calculation is complete */ 02216 if( strcmp(chTime,"LAST") == 0 ) 02217 { 02218 double absval , rel; 02219 long int n; 02220 /* the lines we will find, 02221 * for a sample list of PDR lines look at LineList_PDR_H2.dat 02222 * in the cloudy data dir */ 02223 /* the number of H2 lines */ 02224 const int NLINE_H2 = 31; 02225 /* the number of lines which are not H2 */ 02226 const int NLINE_NOTH_H2 = 5; 02227 /* the labels and wavelengths for the lines that are not H2 */ 02228 char chLabel[NLINE_NOTH_H2][5]= 02229 {"C 2", "O 1", "O 1","C 1", "C 1" }; 02230 double Wl[NLINE_NOTH_H2]= 02231 {157.6 , 63.17 , 145.5 ,609.2 , 369.7 , }; 02232 /* these are wavelengths in microns, conv to Angstroms before call */ 02233 /* >>chng 05 sep 06, many of following wavelengths updated to agree 02234 * with output - apparently not updated when energies changed */ 02235 double Wl_H2[NLINE_H2]= 02236 {2.121 , 02237 28.21 , 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 , 02238 4.408 , 4.180 , 3.996 , 3.845 , 3.724 , 3.626 , 3.547 , 3.485 , 3.437 , 3.403 , 02239 3.381 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603}; 02240 /* print a header for the lines */ 02241 for( n=0; n<NLINE_NOTH_H2; ++n ) 02242 { 02243 fprintf(punch.ipPnunit[ipPun], "%s\t%.2f",chLabel[n] , Wl[n]); 02244 /* get the line, non positive return says didn't find it */ 02245 /* arguments are 4-char label, wavelength, return log total intensity, linear rel inten */ 02246 if( cdLine( chLabel[n] , (realnum)(Wl[n]*1e4) , &absval , &rel ) <= 0 ) 02247 { 02248 fprintf(punch.ipPnunit[ipPun], " did not find\n"); 02249 } 02250 else 02251 { 02252 fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval); 02253 } 02254 } 02255 fprintf(punch.ipPnunit[ipPun], "\n\n\n"); 02256 02257 /* only print the H2 lines if the big molecule is turned on */ 02258 if( h2.lgH2ON ) 02259 { 02260 fprintf(punch.ipPnunit[ipPun], 02261 "Here are some of the H2 Intensities, The first one is the\n" 02262 "1-0 S(0) line and the following ones are the 0-0 S(X)\n" 02263 "lines where X goes from 0 to 29\n\n"); 02264 for( n=0; n<NLINE_H2; ++n ) 02265 { 02266 fprintf(punch.ipPnunit[ipPun], "%s\t%.3f","H2 " , Wl_H2[n]); 02267 /* get the line, non positive return says didn't find it */ 02268 if( cdLine( "H2 " , (realnum)(Wl_H2[n]*1e4) , &absval , &rel ) <= 0 ) 02269 { 02270 fprintf(punch.ipPnunit[ipPun], " did not find\n"); 02271 } 02272 else 02273 { 02274 fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval); 02275 } 02276 } 02277 } 02278 } 02279 } 02280 02281 else if( strcmp( punch.chPunch[ipPun],"LEIS") == 0) 02282 { 02283 if( strcmp(chTime,"LAST") != 0 ) 02284 { 02285 /* get some column densities we shall need */ 02286 double col_ci , col_oi , col_cii, col_heii; 02287 if( cdColm("carb" , 1 , &col_ci ) ) 02288 TotalInsanity(); 02289 if( cdColm("carb" , 2 , &col_cii ) ) 02290 TotalInsanity(); 02291 if( cdColm("oxyg" , 1 , &col_oi ) ) 02292 TotalInsanity(); 02293 if( cdColm("heli" , 2 , &col_heii ) ) 02294 TotalInsanity(); 02295 /* punch Leiden structure - some numbers for the Leiden PDR model comparisons */ 02296 fprintf( punch.ipPnunit[ipPun], 02297 "%.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" 02298 "%.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" 02299 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t" 02300 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t" 02301 "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 02302 /* depth of this point */ 02303 radius.depth_mid_zone, 02304 /* A_V for an extended source */ 02305 0.00, 02306 /* A_V for a point source */ 02307 rfield.extin_mag_V_point, 02308 /* temperature */ 02309 phycon.te , 02310 dense.xIonDense[ipHYDROGEN][0], 02311 hmi.H2_total, 02312 dense.xIonDense[ipCARBON][0], 02313 dense.xIonDense[ipCARBON][1], 02314 dense.xIonDense[ipOXYGEN][0], 02315 findspecies("CO")->hevmol, 02316 findspecies("O2")->hevmol, 02317 findspecies("CH")->hevmol, 02318 findspecies("OH")->hevmol, 02319 dense.eden, 02320 dense.xIonDense[ipHELIUM][1], 02321 dense.xIonDense[ipHYDROGEN][1], 02322 hmi.Hmolec[ipMH3p], 02323 colden.colden[ipCOL_H0], 02324 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s], 02325 col_ci, 02326 col_cii, 02327 col_oi, 02328 findspecies("CO")->hevcol, 02329 findspecies("O2")->hevcol, 02330 findspecies("CH")->hevcol, 02331 findspecies("OH")->hevcol, 02332 colden.colden[ipCOL_elec], 02333 col_heii, 02334 colden.colden[ipCOL_Hp], 02335 colden.colden[ipCOL_H3p], 02336 hmi.H2_Solomon_dissoc_rate_used_H2g , 02337 gv.rate_h2_form_grains_used_total, 02338 hmi.UV_Cont_rel2_Draine_DB96_depth, 02339 /* CO and C dissociation rate */ 02340 CO_findrk("PHOTON,CO=>C,O"), 02341 /* total CI ionization rate */ 02342 ionbal.PhotoRate_Shell[ipCARBON][0][2][0], 02343 /* total heating, erg cm-3 s-1 */ 02344 thermal.htot, 02345 /* total cooling, erg cm-3 s-1 */ 02346 thermal.ctot, 02347 /* GrnP grain photo heating */ 02348 thermal.heating[0][13], 02349 /* grain collisional cooling */ 02350 MAX2(0.,gv.GasCoolColl), 02351 /* grain collisional heating */ 02352 -1.*MIN2(0.,gv.GasCoolColl), 02353 /* COds - CO dissociation heating */ 02354 thermal.heating[0][9], 02355 /* H2dH-Heating due to H2 dissociation */ 02356 hmi.HeatH2Dish_used, 02357 /* H2vH-Heating due to collisions with H2 */ 02358 hmi.HeatH2Dexc_used , 02359 /* ChaT - charge transfer heating */ 02360 thermal.heating[0][24] , 02361 /* cosmic ray heating */ 02362 thermal.heating[1][6] , 02363 /* heating due to atoms of various heavy elements */ 02364 thermal.heating[ipMAGNESIUM][0], 02365 thermal.heating[ipSULPHUR][0], 02366 thermal.heating[ipSILICON][0], 02367 thermal.heating[ipIRON][0], 02368 thermal.heating[ipSODIUM][0], 02369 thermal.heating[ipALUMINIUM][0], 02370 thermal.heating[ipCARBON][0], 02371 TauLines[ipT610].Coll.cool, 02372 TauLines[ipT370].Coll.cool, 02373 TauLines[ipT157].Coll.cool, 02374 TauLines[ipT63].Coll.cool, 02375 TauLines[ipT146].Coll.cool ); 02376 } 02377 } 02378 02379 else if( strcmp( punch.chPunch[ipPun],"LLST") == 0) 02380 { 02381 /* punch linelist command - do on last iteration */ 02382 if( strcmp(chTime,"LAST") == 0 ) 02383 { 02384 fprintf( punch.ipPnunit[ipPun], 02385 "%li" , iteration ); 02386 02387 /* -1 is flag saying that this punch command was not set */ 02388 if( punch.nLineList[ipPun] < 0 ) 02389 TotalInsanity(); 02390 02391 /* loop over all lines in the file we read */ 02392 for( j=0; j<punch.nLineList[ipPun]; ++j ) 02393 { 02394 double relative , absolute, PrtQuantity; 02395 if( (cdLine( punch.chLineListLabel[ipPun][j] , 02396 punch.wlLineList[ipPun][j] , 02397 &relative , &absolute ))<=0 ) 02398 { 02399 if( !h2.lgH2ON && strncmp( punch.chLineListLabel[ipPun][j] , "H2 " , 4 )==0 ) 02400 { 02401 static bool lgMustPrintFirstTime = true; 02402 if( lgMustPrintFirstTime ) 02403 { 02404 /* it's an H2 line and H2 is not being done - ignore it */ 02405 fprintf( ioQQQ,"Did not find an H2 line, the large model is not " 02406 "included, so I will ignore it. Log intensity set to -30.\n" ); 02407 fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n"); 02408 lgMustPrintFirstTime = false; 02409 } 02410 relative = -30.f; 02411 absolute = -30.f; 02412 } 02413 else if( lgAbort ) 02414 { 02415 /* we are in abort mode */ 02416 relative = -30.f; 02417 absolute = -30.f; 02418 } 02419 else 02420 { 02421 fprintf(ioQQQ,"DISASTER - did not find a line in the Line List table\n"); 02422 cdEXIT(EXIT_FAILURE); 02423 } 02424 } 02425 02426 /* options to do either relative or absolute intensity 02427 * default is relative, is absolute keyword on line then 02428 * punarg set to 1 */ 02429 /* straight line intensities */ 02430 if( punch.punarg[ipPun][0] > 0 ) 02431 PrtQuantity = absolute; 02432 else 02433 PrtQuantity = relative; 02434 /* if taking ratio print every other line as ratio 02435 * with previous line */ 02436 if( punch.lgLineListRatio[ipPun] ) 02437 { 02438 /* do line pair ratios */ 02439 static double SaveQuantity = 0; 02440 if( is_odd(j) ) 02441 fprintf( punch.ipPnunit[ipPun], "\t%.3e" , 02442 SaveQuantity / SDIV( PrtQuantity ) ); 02443 else 02444 SaveQuantity = PrtQuantity; 02445 } 02446 else 02447 { 02448 fprintf( punch.ipPnunit[ipPun], "\t%.3e" , PrtQuantity ); 02449 } 02450 } 02451 fprintf( punch.ipPnunit[ipPun], "\n" ); 02452 } 02453 } 02454 02455 else if( strcmp( punch.chPunch[ipPun],"RCO ") == 0) 02456 { 02457 /* punch chemistry CO rates command */ 02458 if( strcmp(chTime,"LAST") != 0 ) 02459 { 02460 CO_punch_mol(punch.ipPnunit[ipPun],"CO",NULL,radius.depth_mid_zone); 02461 } 02462 } 02463 02464 /*>>chng 06 jul 13, Nick Abel add OH rates */ 02465 else if( strcmp( punch.chPunch[ipPun],"ROH ") == 0) 02466 { 02467 /* punch chemistry OH rates command */ 02468 if( strcmp(chTime,"LAST") != 0 ) 02469 { 02470 CO_punch_mol(punch.ipPnunit[ipPun],"OH",NULL,radius.depth_mid_zone); 02471 } 02472 } 02473 else if( strcmp(punch.chPunch[ipPun],"MAP ") == 0 ) 02474 { 02475 /* do the map now if we are at the zone, or if this 02476 * is the LAST call to this routine and map not done yet */ 02477 if( !hcmap.lgMapDone && 02478 (nzone == hcmap.MapZone || strcmp(chTime,"LAST") == 0 ) ) 02479 { 02480 lgTlkSav = called.lgTalk; 02481 called.lgTalk = true; 02482 hcmap.lgMapBeingDone = true; 02483 map_do(punch.ipPnunit[ipPun] , " map"); 02484 called.lgTalk = lgTlkSav; 02485 } 02486 } 02487 02488 else if( strcmp(punch.chPunch[ipPun],"MOLE") == 0 ) 02489 { 02490 if( strcmp(chTime,"LAST") != 0 ) 02491 { 02492 static bool lgMustPrintHeader = true; 02493 if( lgMustPrintHeader ) 02494 { 02495 lgMustPrintHeader = false; 02496 fprintf( punch.ipPnunit[ipPun], 02497 "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate"); 02498 02499 for(i=0; i<N_H_MOLEC; ++i ) 02500 { 02501 fprintf( punch.ipPnunit[ipPun], "\t%s", hmi.chLab[i] ); 02502 } 02503 for(i=0; i<mole.num_comole_calc; ++i ) 02504 { 02505 fprintf( punch.ipPnunit[ipPun], "\t%s",COmole[i]->label ); 02506 } 02507 fprintf( punch.ipPnunit[ipPun], "\n"); 02508 } 02509 /* molecules, especially for PDR, first give radius */ 02510 fprintf( punch.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone ); 02511 02512 /* visual extinction of point source (star)*/ 02513 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point); 02514 02515 /* visual extinction of an extended source (like a PDR)*/ 02516 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended); 02517 02518 /* carbon monoxide photodissociation rate */ 02519 fprintf( punch.ipPnunit[ipPun], "%.5e\t" , CO_findrk("PHOTON,CO=>C,O") ); 02520 02521 /* carbon recombination rate */ 02522 fprintf( punch.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] ); 02523 02524 /* now do all the H molecules */ 02525 for(j=0; j<N_H_MOLEC; ++j ) 02526 { 02527 fprintf(punch.ipPnunit[ipPun],"\t%.2e",hmi.Hmolec[j] ); 02528 } 02529 02530 /* now all the C molecules */ 02531 for(j=0; j<mole.num_comole_calc; ++j ) 02532 { 02533 fprintf(punch.ipPnunit[ipPun],"\t%.2e",COmole[j]->hevmol ); 02534 } 02535 fprintf(punch.ipPnunit[ipPun],"\n"); 02536 } 02537 } 02538 02539 else if( strcmp(punch.chPunch[ipPun],"OPAC") == 0 ) 02540 { 02541 /* punch opacity - routine will parse which type of opacity punch to do */ 02542 if( strcmp(chTime,"LAST") == 0 ) 02543 punch_opacity(punch.ipPnunit[ipPun],ipPun); 02544 } 02545 02546 /* punch coarse optical depths command */ 02547 else if( strcmp(punch.chPunch[ipPun],"OPTc") == 0 ) 02548 { 02549 if( strcmp(chTime,"LAST") == 0 ) 02550 { 02551 for( j=0; j < rfield.nflux; j++ ) 02552 { 02553 fprintf( punch.ipPnunit[ipPun], 02554 "%12.4e\t%.3e\t%12.4e\t%.3e\n", 02555 AnuUnit(rfield.AnuOrg[j]), 02556 opac.TauAbsFace[j]+opac.TauScatFace[j], 02557 opac.TauAbsFace[j], 02558 opac.TauScatFace[j] ); 02559 } 02560 } 02561 } 02562 02563 /* punch fine optical depths command */ 02564 else if( strcmp(punch.chPunch[ipPun],"OPTf") == 0 ) 02565 { 02566 if( strcmp(chTime,"LAST") == 0 ) 02567 { 02568 long nu_hi , nskip; 02569 if( punch.punarg[ipPun][0] > 0. ) 02570 /* possible lower bounds to energy range - will be zero if not set */ 02571 j = ipFineCont( punch.punarg[ipPun][0] ); 02572 else 02573 j = 0; 02574 02575 /* upper limit */ 02576 if( punch.punarg[ipPun][1]> 0. ) 02577 nu_hi = ipFineCont( punch.punarg[ipPun][1]); 02578 else 02579 nu_hi = rfield.nfine; 02580 02581 /* we will bring nskip cells together into one printed 02582 * number to make output smaller - default is 10 */ 02583 nskip = (long)punch.punarg[ipPun][2]; 02584 do 02585 { 02586 realnum sum1 = rfield.fine_opt_depth[j]; 02587 realnum sum2 = rfield.fine_opac_zone[j]; 02588 /* want to report the central wavelength of the cell */ 02589 realnum xnu = rfield.fine_anu[j]; 02590 for( jj=1; jj<nskip; ++jj ) 02591 { 02592 sum1 += rfield.fine_opt_depth[j+jj]; 02593 sum2 += rfield.fine_opac_zone[j+jj]; 02594 xnu += rfield.fine_anu[j+jj]; 02595 } 02596 if( sum2 > 0. ) 02597 fprintf( punch.ipPnunit[ipPun], 02598 "%12.6e\t%.3e\t%.3e\n", 02599 AnuUnit(xnu/nskip), 02600 sum1/nskip , 02601 sum2/nskip); 02602 j += nskip; 02603 }while( j < nu_hi ); 02604 } 02605 } 02606 02607 else if( strcmp(punch.chPunch[ipPun]," OTS") == 0 ) 02608 { 02609 ConMax = 0.; 02610 xLinMax = 0.; 02611 opConSum = 0.; 02612 opLinSum = 0.; 02613 ipLinMax = 1; 02614 ipConMax = 1; 02615 02616 for( j=0; j < rfield.nflux; j++ ) 02617 { 02618 opConSum += rfield.otscon[j]*opac.opacity_abs[j]; 02619 opLinSum += rfield.otslin[j]*opac.opacity_abs[j]; 02620 if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax ) 02621 { 02622 xLinMax = rfield.otslin[j]*opac.opacity_abs[j]; 02623 ipLinMax = j+1; 02624 } 02625 if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax ) 02626 { 02627 ConMax = rfield.otscon[j]*opac.opacity_abs[j]; 02628 ipConMax = j+1; 02629 } 02630 } 02631 fprintf( punch.ipPnunit[ipPun], 02632 "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n", 02633 opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1] 02634 , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1] 02635 , rfield.anu[ipConMax-1], ConMax ); 02636 } 02637 02638 else if( strcmp(punch.chPunch[ipPun],"OVER") == 0 ) 02639 { 02640 /* punch overview 02641 * this is the floor for the smallest ionization fractions printed */ 02642 double toosmall = SMALLFLOAT , 02643 hold; 02644 02645 /* overview of model results, 02646 * depth, te, hden, eden, ion fractions H, He, c, O */ 02647 if( strcmp(chTime,"LAST") != 0 ) 02648 { 02649 02650 /* print the depth */ 02651 fprintf( punch.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone ); 02652 02653 /* temperature, heating */ 02654 if(dynamics.Cool > dynamics.Heat) 02655 { 02656 fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f", 02657 log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat)) ); 02658 } 02659 else 02660 { 02661 double diff = fabs(thermal.htot-dynamics.Cool); 02662 fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f", 02663 log10(phycon.te), log10(SDIV(diff)) ); 02664 } 02665 02666 /* hydrogen and electron densities */ 02667 fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f", 02668 log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) ); 02669 02670 /* molecular fraction of hydrogen */ 02671 fprintf( punch.ipPnunit[ipPun], "\t%.4f", 02672 /*log10(MAX2(toosmall,2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN])) );*/ 02673 log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) ); 02674 02675 /* ionization fractions of hydrogen */ 02676 fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f", 02677 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])), 02678 log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) ); 02679 02680 /* ionization fractions of helium */ 02681 for( j=1; j <= 3; j++ ) 02682 { 02683 double arg1 = SDIV(dense.gas_phase[ipHELIUM]); 02684 arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 ); 02685 fprintf( punch.ipPnunit[ipPun], "\t%.4f", 02686 log10(arg1) ); 02687 } 02688 02689 /* carbon monoxide molecular fraction of CO */ 02690 hold = SDIV(dense.gas_phase[ipCARBON]); 02691 hold = findspecies("CO")->hevmol/hold; 02692 hold = MAX2(toosmall, hold ); 02693 fprintf( punch.ipPnunit[ipPun], "\t%.4f", log10(hold) ); 02694 02695 /* ionization fractions of carbon */ 02696 for( j=1; j <= 4; j++ ) 02697 { 02698 hold = SDIV(dense.gas_phase[ipCARBON]); 02699 hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold); 02700 fprintf( punch.ipPnunit[ipPun], "\t%.4f", 02701 log10(hold) ); 02702 } 02703 02704 /* ionization fractions of oxygen */ 02705 for( j=1; j <= 6; j++ ) 02706 { 02707 hold = SDIV(dense.gas_phase[ipOXYGEN]); 02708 hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold); 02709 fprintf( punch.ipPnunit[ipPun], "\t%.4f", 02710 log10(hold) ); 02711 } 02712 02713 /* visual extinction of point source (star)*/ 02714 fprintf( punch.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point); 02715 02716 /* visual extinction of an extended source (like a PDR)*/ 02717 fprintf( punch.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended); 02718 } 02719 } 02720 02721 else if( strcmp(punch.chPunch[ipPun]," PDR") == 0 ) 02722 { 02723 /* this is the punch PDR command */ 02724 if( strcmp(chTime,"LAST") != 0 ) 02725 { 02726 /* convert optical depth at wavelength of V filter 02727 * into magnitudes of extinction */ 02728 /* >>chyng 03 feb 25, report extinction to illuminated face, 02729 * rather than total extinction which included far side when 02730 * sphere was set */ 02731 /*av = opac.TauTotalGeo[0][rfield.ipV_filter-1]*1.08574;*/ 02732 02733 fprintf( punch.ipPnunit[ipPun], 02734 "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t", 02735 radius.depth_mid_zone, 02736 /* total hydrogen column density, all forms */ 02737 colden.colden[ipCOL_HTOT], 02738 phycon.te, 02739 /* fraction of H that is atomic */ 02740 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN], 02741 /* ratio of n(H2) to total H, == 0.5 when fully molecular */ 02742 2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN], 02743 2.*hmi.Hmolec[ipMH2s]/dense.gas_phase[ipHYDROGEN], 02744 /* atomic to total carbon */ 02745 dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]), 02746 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 02747 findspecies("H2O")->hevmol/SDIV(dense.gas_phase[ipOXYGEN]), 02748 /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */ 02749 hmi.UV_Cont_rel2_Habing_TH85_depth); 02750 02751 /* visual extinction due to dust alone, of point source (star)*/ 02752 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point); 02753 02754 /* visual extinction due to dust alone, of an extended source (like a PDR)*/ 02755 fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended); 02756 02757 /* visual extinction (all sources) of a point source (like a PDR)*/ 02758 fprintf( punch.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] ); 02759 } 02760 } 02761 02762 else if( strcmp(punch.chPunch[ipPun],"PHYS") == 0 ) 02763 { 02764 if( strcmp(chTime,"LAST") != 0 ) 02765 { 02766 /* punch physical conditions */ 02767 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 02768 radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], 02769 dense.eden, thermal.htot, wind.AccelTot, geometry.FillFac ); 02770 } 02771 } 02772 02773 else if( strcmp(punch.chPunch[ipPun],"PRES") == 0 ) 02774 { 02775 /* the punch pressure command */ 02776 if( strcmp(chTime,"LAST") != 0 ) 02777 { 02778 fprintf( punch.ipPnunit[ipPun], 02779 "%.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", 02780 /*A 1 #P depth */ 02781 radius.depth_mid_zone, 02782 /*B 2 Pcorrect */ 02783 pressure.PresTotlCorrect, 02784 /*C 3 Pcurrent */ 02785 pressure.PresTotlCurr, 02786 /*D 4 Pln + pintg 02787 * >>chng 06 apr 19, subtract pinzon the acceleration added in this zone 02788 * since is not total at outer edge of zone, above is at inner edge */ 02789 pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon, 02790 /*E 5 pgas (0) */ 02791 pressure.PresTotlInit, 02792 /*F 6 Pgas */ 02793 pressure.PresGasCurr, 02794 /*G 7 Pram */ 02795 pressure.PresRamCurr, 02796 /*H 8 P rad in lines */ 02797 pressure.pres_radiation_lines_curr, 02798 /*I 9 Pinteg subtract continuum rad pres which has already been added on */ 02799 pressure.PresInteg - pressure.pinzon, 02800 /*J 10 V(wind km/s) wind speed in km/s */ 02801 wind.windv/1e5, 02802 /*K cad(km/s) sound speed in km/s */ 02803 timesc.sound_speed_adiabatic/1e5, 02804 /* the magnetic pressure */ 02805 magnetic.pressure , 02806 /* the local turbulent velocity in km/s */ 02807 DoppVel.TurbVel/1e5 , 02808 /* turbulent pressure */ 02809 pressure.PresTurbCurr*DoppVel.lgTurb_pressure , 02810 TorF(conv.lgConvPres) ); 02811 } 02812 } 02813 02814 else if( punch.chPunch[ipPun][0]=='R' ) 02815 { 02816 /* work around internal limits to Microsoft vs compiler */ 02817 if( strcmp(punch.chPunch[ipPun],"RADI") == 0 ) 02818 { 02819 /* punch radius information for all zones */ 02820 if( strcmp(chTime,"LAST") != 0 ) 02821 { 02822 fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n", 02823 nzone, radius.Radius_mid_zone, radius.depth_mid_zone, 02824 radius.drad ); 02825 } 02826 } 02827 02828 else if( strcmp(punch.chPunch[ipPun],"RADO") == 0 ) 02829 { 02830 /* punch radius information for only the last zone */ 02831 if( strcmp(chTime,"LAST") == 0 ) 02832 { 02833 fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n", 02834 nzone, radius.Radius_mid_zone, radius.depth_mid_zone, 02835 radius.drad ); 02836 } 02837 } 02838 02839 else if( strcmp(punch.chPunch[ipPun],"RESU") == 0 ) 02840 { 02841 /* punch results of the calculation */ 02842 if( strcmp(chTime,"LAST") == 0 ) 02843 punResults(punch.ipPnunit[ipPun]); 02844 } 02845 else 02846 { 02847 /* this can't happen */ 02848 TotalInsanity(); 02849 } 02850 } 02851 02852 else if( strcmp(punch.chPunch[ipPun],"SECO") == 0 ) 02853 { 02854 /* punch secondary ionization */ 02855 if( strcmp(chTime,"LAST") != 0 ) 02856 fprintf(punch.ipPnunit[ipPun], 02857 "%.5e\t%.3e\t%.3e\t%.3e\n", 02858 radius.depth , 02859 secondaries.csupra[ipHYDROGEN][0], 02860 secondaries.csupra[ipHYDROGEN][0]*2.02, 02861 secondaries.x12tot ); 02862 } 02863 02864 else if( strcmp(punch.chPunch[ipPun],"SOUS") == 0 ) 02865 { 02866 /* full spectrum of continuum source function at 1 depth 02867 * command was "punch source spectrum" */ 02868 if( strcmp(chTime,"LAST") != 0 ) 02869 { 02870 limit = MIN2(rfield.ipMaxBolt,rfield.nflux); 02871 for( j=0; j < limit; j++ ) 02872 { 02873 fprintf( punch.ipPnunit[ipPun], 02874 "%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 02875 rfield.anu[j], 02876 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j], 02877 opac.opacity_abs[j], 02878 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j]), 02879 rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j])/plankf(j) ); 02880 } 02881 } 02882 } 02883 02884 else if( strcmp(punch.chPunch[ipPun],"SOUD") == 0 ) 02885 { 02886 /* parts of continuum source function vs depth 02887 * command was punch source function depth */ 02888 j = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] + 2; 02889 fprintf( punch.ipPnunit[ipPun], 02890 "%.4e\t%.4e\t%.4e\t%.4e\n", 02891 opac.TauAbsFace[j-1], 02892 rfield.ConEmitLocal[nzone][j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]), 02893 rfield.otscon[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1], 02894 rfield.otscon[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/opac.opacity_abs[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] ); 02895 } 02896 02897 /* this is punch special option */ 02898 else if( strcmp(punch.chPunch[ipPun],"SPEC") == 0 ) 02899 { 02900 PunchSpecial(punch.ipPnunit[ipPun],chTime); 02901 } 02902 02903 else if( strcmp(punch.chPunch[ipPun],"TEGR") == 0 ) 02904 { 02905 /* punch history of heating and cooling */ 02906 if( strcmp(chTime,"LAST") != 0 ) 02907 { 02908 fprintf( punch.ipPnunit[ipPun], " " ); 02909 for( j=0; j < NGRID; j++ ) 02910 { 02911 fprintf( punch.ipPnunit[ipPun], 02912 "%4ld%11.3e%11.3e%11.3e\n", 02913 thermal.nZonGrid[j], 02914 thermal.TeGrid[j], 02915 thermal.HtGrid[j], 02916 thermal.ClGrid[j] ); 02917 } 02918 } 02919 } 02920 02921 else if( strcmp(punch.chPunch[ipPun],"TEMP") == 0 ) 02922 { 02923 static double deriv_old=-1; 02924 double deriv=-1. , deriv_sec; 02925 /* temperature and its derivatives */ 02926 fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e", 02927 radius.depth_mid_zone, 02928 phycon.te, 02929 thermal.dCooldT ); 02930 /* if second zone then have one deriv */ 02931 if( nzone >1 ) 02932 { 02933 deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad; 02934 fprintf( punch.ipPnunit[ipPun], "\t%.2e", deriv ); 02935 /* if third zone then have second deriv */ 02936 if( nzone > 2 ) 02937 { 02938 deriv_sec = (deriv-deriv_old)/ radius.drad; 02939 fprintf( punch.ipPnunit[ipPun], "\t%.2e", 02940 deriv_sec ); 02941 } 02942 deriv_old = deriv; 02943 } 02944 fprintf( punch.ipPnunit[ipPun], "\n"); 02945 } 02946 02947 /* time dependent model */ 02948 else if( strcmp(punch.chPunch[ipPun],"TIMD") == 0 ) 02949 { 02950 if( strcmp(chTime,"LAST") == 0 ) 02951 DynaPunchTimeDep( punch.ipPnunit[ipPun] , "END" ); 02952 } 02953 02954 /* execution time per zone */ 02955 else if( strcmp(punch.chPunch[ipPun],"XTIM") == 0 ) 02956 { 02957 static double ElapsedTime , ZoneTime; 02958 if( nzone<=1 ) 02959 { 02960 ElapsedTime = cdExecTime(); 02961 ZoneTime = 0.; 02962 } 02963 else 02964 { 02965 double t = cdExecTime(); 02966 ZoneTime = t - ElapsedTime; 02967 ElapsedTime = t; 02968 } 02969 02970 /* zone, time for this zone, elapsed time */ 02971 fprintf( punch.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n", 02972 nzone, ZoneTime , ElapsedTime ); 02973 } 02974 02975 else if( strcmp(punch.chPunch[ipPun],"TPRE") == 0 ) 02976 { 02977 /* temperature and its predictors, turned on with punch tprid */ 02978 fprintf( punch.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n", 02979 nzone, phycon.TeInit, phycon.TeProp, phycon.te, 02980 (phycon.TeProp- phycon.te)/phycon.te ); 02981 } 02982 02983 else if( strcmp(punch.chPunch[ipPun],"WIND") == 0 ) 02984 { 02985 /* wind velocity, radiative acceleration, and ratio total 02986 * to electron scattering acceleration */ 02987 /* first test only punch last zone */ 02988 if( (punch.punarg[ipPun][0] == 0 && strcmp(chTime,"LAST") == 0) 02989 || 02990 /* this test punch all zones */ 02991 (punch.punarg[ipPun][0] == 1 && strcmp(chTime,"LAST") != 0 ) ) 02992 { 02993 fprintf( punch.ipPnunit[ipPun], 02994 "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 02995 radius.Radius_mid_zone, 02996 radius.depth_mid_zone, 02997 wind.windv, 02998 wind.AccelTot, 02999 wind.AccelLine, 03000 wind.AccelCont , 03001 wind.fmul , 03002 wind.AccelGravity ); 03003 } 03004 } 03005 03006 else if( strcmp(punch.chPunch[ipPun],"XATT") == 0 ) 03007 { 03008 /* attenuated incident continuum */ 03009 ASSERT( grid.lgOutputTypeOn[2] ); 03010 03011 if( strcmp(chTime,"LAST") == 0 ) 03012 { 03013 if( grid.lgGrid ) 03014 punchFITSfile( punch.ipPnunit[ipPun], 2 ); 03015 else 03016 { 03017 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03018 cdEXIT(EXIT_FAILURE); 03019 } 03020 } 03021 } 03022 else if( strcmp(punch.chPunch[ipPun],"XRFI") == 0 ) 03023 { 03024 /* reflected incident continuum */ 03025 ASSERT( grid.lgOutputTypeOn[3] ); 03026 03027 if( strcmp(chTime,"LAST") == 0 ) 03028 { 03029 if( grid.lgGrid ) 03030 punchFITSfile( punch.ipPnunit[ipPun], 3 ); 03031 else 03032 { 03033 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03034 cdEXIT(EXIT_FAILURE); 03035 } 03036 } 03037 } 03038 else if( strcmp(punch.chPunch[ipPun],"XINC") == 0 ) 03039 { 03040 /* incident continuum */ 03041 ASSERT( grid.lgOutputTypeOn[1] ); 03042 03043 if( strcmp(chTime,"LAST") == 0 ) 03044 { 03045 if( grid.lgGrid ) 03046 punchFITSfile( punch.ipPnunit[ipPun], 1 ); 03047 else 03048 { 03049 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03050 cdEXIT(EXIT_FAILURE); 03051 } 03052 } 03053 } 03054 else if( strcmp(punch.chPunch[ipPun],"XDFR") == 0 ) 03055 { 03056 /* reflected diffuse continuous emission */ 03057 ASSERT( grid.lgOutputTypeOn[5] ); 03058 03059 if( strcmp(chTime,"LAST") == 0 ) 03060 { 03061 if( grid.lgGrid ) 03062 punchFITSfile( punch.ipPnunit[ipPun], 5 ); 03063 else 03064 { 03065 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03066 cdEXIT(EXIT_FAILURE); 03067 } 03068 } 03069 } 03070 else if( strcmp(punch.chPunch[ipPun],"XDFO") == 0 ) 03071 { 03072 /* diffuse continuous emission outward */ 03073 ASSERT( grid.lgOutputTypeOn[4] ); 03074 03075 if( strcmp(chTime,"LAST") == 0 ) 03076 { 03077 if( grid.lgGrid ) 03078 punchFITSfile( punch.ipPnunit[ipPun], 4 ); 03079 else 03080 { 03081 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03082 cdEXIT(EXIT_FAILURE); 03083 } 03084 } 03085 } 03086 else if( strcmp(punch.chPunch[ipPun],"XLNR") == 0 ) 03087 { 03088 /* reflected lines */ 03089 ASSERT( grid.lgOutputTypeOn[7] ); 03090 03091 if( strcmp(chTime,"LAST") == 0 ) 03092 { 03093 if( grid.lgGrid ) 03094 punchFITSfile( punch.ipPnunit[ipPun], 7 ); 03095 else 03096 { 03097 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03098 cdEXIT(EXIT_FAILURE); 03099 } 03100 } 03101 } 03102 else if( strcmp(punch.chPunch[ipPun],"XLNO") == 0 ) 03103 { 03104 /* outward lines */ 03105 ASSERT( grid.lgOutputTypeOn[6] ); 03106 03107 if( strcmp(chTime,"LAST") == 0 ) 03108 { 03109 if( grid.lgGrid ) 03110 punchFITSfile( punch.ipPnunit[ipPun], 6 ); 03111 else 03112 { 03113 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03114 cdEXIT(EXIT_FAILURE); 03115 } 03116 } 03117 } 03118 else if( strcmp(punch.chPunch[ipPun],"XREF") == 0 ) 03119 { 03120 /* total reflected, lines and continuum */ 03121 ASSERT( grid.lgOutputTypeOn[9] ); 03122 03123 if( strcmp(chTime,"LAST") == 0 ) 03124 { 03125 if( grid.lgGrid ) 03126 punchFITSfile( punch.ipPnunit[ipPun], 9 ); 03127 else 03128 { 03129 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03130 cdEXIT(EXIT_FAILURE); 03131 } 03132 } 03133 } 03134 else if( strcmp(punch.chPunch[ipPun],"XTRN") == 0 ) 03135 { 03136 /* total outward, lines and continuum */ 03137 ASSERT( grid.lgOutputTypeOn[8] ); 03138 03139 if( strcmp(chTime,"LAST") == 0 ) 03140 { 03141 if( grid.lgGrid ) 03142 punchFITSfile( punch.ipPnunit[ipPun], 8 ); 03143 else 03144 { 03145 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03146 cdEXIT(EXIT_FAILURE); 03147 } 03148 } 03149 } 03150 else if( strcmp(punch.chPunch[ipPun],"XSPM") == 0 ) 03151 { 03152 /* exp(-tau) to the illuminated face */ 03153 ASSERT( grid.lgOutputTypeOn[10] ); 03154 03155 if( strcmp(chTime,"LAST") == 0 ) 03156 { 03157 if( grid.lgGrid ) 03158 punchFITSfile( punch.ipPnunit[ipPun], 10 ); 03159 else 03160 { 03161 fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" ); 03162 cdEXIT(EXIT_FAILURE); 03163 } 03164 } 03165 } 03166 /* there are a few "punch" commands that are handled elsewhere 03167 * punch dr is an example. These will have lgReadPunch set 03168 * false */ 03169 else if( punch.lgRealPunch[ipPun] ) 03170 { 03171 /* this is insanity, internal flag set in ParsePunch not seen here */ 03172 fprintf( ioQQQ, " PROBLEM DISASTER PunchDo does not recognize flag %4.4s set by ParsePunch. This is impossible.\n", 03173 punch.chPunch[ipPun] ); 03174 TotalInsanity(); 03175 } 03176 03177 /* print special hash string to separate out various iterations 03178 * chTime is LAST on last iteration 03179 * punch.lgHashEndIter flag is true by default, set false 03180 * with "no hash" keyword on punch command 03181 * punch.lg_separate_iterations is true by default, set false 03182 * when punch time dependent calcs since do not want special 03183 * character between time steps 03184 * grid.lgGrid is only true when doing a grid of calculations */ 03185 if( strcmp(chTime,"LAST") == 0 && 03186 !(iterations.lgLastIt && !grid.lgGrid ) && 03187 punch.lgHashEndIter[ipPun] && 03188 punch.lg_separate_iterations[ipPun] && 03189 !punch.lgFITS[ipPun] ) 03190 { 03191 if( dynamics.lgStatic && strcmp( punch.chHashString , "TIME_DEP" )==0 ) 03192 { 03193 fprintf( punch.ipPnunit[ipPun], "\"time=%f\n", 03194 dynamics.time_elapsed ); 03195 } 03196 else 03197 { 03198 fprintf( punch.ipPnunit[ipPun], "%s\n", 03199 punch.chHashString ); 03200 } 03201 } 03202 if( punch.lgFLUSH ) 03203 fflush( punch.ipPnunit[ipPun] ); 03204 } 03205 } 03206 return; 03207 } 03208 03209 /*PunLineIntensity produce the 'punch lines intensity' output */ 03210 STATIC void PunLineIntensity(FILE * ioPUN) 03211 { 03212 char chCard[INPUT_LINE_LENGTH]; 03213 bool lgEOF; 03214 long int i; 03215 03216 DEBUG_ENTRY( "PunLineIntensity()" ); 03217 03218 /* used to punch out all the emission line intensities 03219 * first initialize the line image reader */ 03220 03221 /* start the file with the input commands */ 03222 input_init(); 03223 fprintf( ioPUN, "**********************************************************************************************************************************\n" ); 03224 03225 lgEOF = false; 03226 while( !lgEOF ) 03227 { 03228 input_readarray(chCard,&lgEOF); 03229 if( !lgEOF ) 03230 { 03231 fprintf( ioPUN, "%s\n", chCard ); 03232 } 03233 } 03234 03235 /* now print any cautions or warnings */ 03236 cdWarnings( ioPUN); 03237 cdCautions( ioPUN); 03238 fprintf( ioPUN, "zone=%5ld\n", nzone ); 03239 fprintf( ioPUN, "**********************************************************************************************************************************\n" ); 03240 fprintf( ioPUN, "begin emission lines\n" ); 03241 03242 /* >>chng 96 jul 03, only punch non-zero intensities */ 03243 PunResults1Line(ioPUN," ",0,0.,"Start"); 03244 for( i=0; i < LineSave.nsum; i++ ) 03245 { 03246 if( LineSv[i].sumlin[LineSave.lgLineEmergent] > 0. ) 03247 { 03248 PunResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength, 03249 LineSv[i].sumlin[LineSave.lgLineEmergent], "Line "); 03250 } 03251 } 03252 03253 PunResults1Line(ioPUN," ",0,0.,"Flush"); 03254 03255 fprintf( ioPUN, " \n" ); 03256 fprintf( ioPUN, "**********************************************************************************************************************************\n" ); 03257 03258 return; 03259 } 03260 03261 /* lgPunchOpticalDepths true says punch optical depths */ 03262 static bool lgPopsFirstCall , lgPunchOpticalDepths; 03263 03264 /*PunchLineStuff punch optical depths or source functions for all transferred lines */ 03265 STATIC void PunchLineStuff( 03266 FILE * ioPUN, 03267 const char *chJob , 03268 realnum xLimit ) 03269 { 03270 03271 long int nelem , ipLo , ipHi , i , ipISO; 03272 long index = 0; 03273 static bool lgFirst=true; 03274 03275 DEBUG_ENTRY( "PunchLineStuff()" ); 03276 03277 /*find out which job this is and set a flag to use later */ 03278 if( strcmp( &*chJob , "optical" ) == 0 ) 03279 { 03280 /* punch line optical depths */ 03281 lgPunchOpticalDepths = true; 03282 lgPopsFirstCall = false; 03283 } 03284 else if( strcmp( &*chJob , "populat" ) == 0 ) 03285 { 03286 lgPunchOpticalDepths = false; 03287 /* level population information */ 03288 if( lgFirst ) 03289 { 03290 lgPopsFirstCall = true; 03291 fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n"); 03292 lgFirst = false; 03293 } 03294 else 03295 { 03296 lgPopsFirstCall = false; 03297 } 03298 } 03299 else 03300 { 03301 fprintf( ioQQQ, " insane job in PunchLineStuff =%s\n", 03302 &*chJob ); 03303 cdEXIT(EXIT_FAILURE); 03304 } 03305 03306 /* loop over all lines, calling put1Line to create info (routine located below) */ 03307 /* hydrogen like lines */ 03308 /* >>chng 02 may 16, had been explicit H and He-like loops */ 03309 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 03310 { 03311 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 03312 { 03313 if( dense.lgElmtOn[nelem] ) 03314 { 03315 /* 06 aug 28, from numLevels_max to _local. */ 03316 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 03317 { 03318 for( ipLo=0; ipLo <ipHi; ipLo++ ) 03319 { 03320 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 03321 continue; 03322 03323 ++index; 03324 pun1Line( &Transitions[ipISO][nelem][ipHi][ipLo] , ioPUN , xLimit , index , dense.xIonDense[nelem][nelem+1-ipISO] ); 03325 } 03326 } 03327 /* also do extra Lyman lines if optical depths are to be done, 03328 * these are line that are included only for absorption, not in the 03329 * model atoms */ 03330 if( lgPunchOpticalDepths ) 03331 { 03332 /* >>chng 02 aug 23, for he-like, had starting on much too high a level since 03333 * index was number of levels - caught by Adrian Turner */ 03334 /* now output extra line lines, starting one above those already done above */ 03335 /*for( ipHi=iso.numLevels_max[ipISO][nelem]; ipHi < iso.nLyman[ipISO]; ipHi++ )*/ 03336 /* 06 aug 28, from numLevels_max to _local. */ 03337 for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1].n+1; ipHi < iso.nLyman[ipISO]; ipHi++ ) 03338 { 03339 ++index; 03340 pun1Line( &ExtraLymanLines[ipISO][nelem][ipHi] , ioPUN , xLimit , index , dense.xIonDense[nelem][nelem+1-ipISO]); 03341 } 03342 } 03343 } 03344 } 03345 } 03346 03347 /* index from 1 to skip over dummy line */ 03348 for( i=1; i < nLevel1; i++ ) 03349 { 03350 ++index; 03351 pun1Line( &TauLines[i] , ioPUN , xLimit , index , 1. ); 03352 } 03353 03354 for( i=0; i < nWindLine; i++ ) 03355 { 03356 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 03357 { 03358 ++index; 03359 pun1Line( &TauLine2[i] , ioPUN , xLimit , index , 1. ); 03360 } 03361 } 03362 03363 for( i=0; i < nUTA; i++ ) 03364 { 03365 ++index; 03366 pun1Line( &UTALines[i] , ioPUN , xLimit , index , 1. ); 03367 } 03368 03369 03370 /* do optical depths of FeII lines, if large atom is turned on */ 03371 FeIIPunchLineStuff( ioPUN , xLimit , index); 03372 03373 /* punch optical depths of H2 lines */ 03374 H2_PunchLineStuff( ioPUN , xLimit , index); 03375 03376 /* C12O16 */ 03377 for( i=0; i < nCORotate; i++ ) 03378 { 03379 ++index; 03380 pun1Line( &C12O16Rotate[i] , ioPUN , xLimit , index , 1. ); 03381 } 03382 03383 /* C13O16 */ 03384 for( i=0; i < nCORotate; i++ ) 03385 { 03386 ++index; 03387 pun1Line( &C13O16Rotate[i] , ioPUN , xLimit , index , 1. ); 03388 } 03389 /*fprintf(ioPUN, "##################################\n"); */ 03390 fprintf( ioPUN , "%s\n",punch.chHashString ); 03391 return; 03392 } 03393 03394 /*pun1Line called by PunchLineStuff to produce output for one line */ 03395 void pun1Line( transition * t , FILE * ioPUN , realnum xLimit , long index , double abundance) 03396 { 03397 03398 if( lgPunchOpticalDepths ) 03399 { 03400 /* optical depths, no special first time, only print them */ 03401 if( t->Emis->TauIn >= xLimit ) 03402 { 03403 /* label like "C 4" or "H 1" */ 03404 fprintf( ioPUN , "%2.2s%2.2s\t", 03405 elementnames.chElementSym[t->Hi->nelem-1] , 03406 elementnames.chIonStage[t->Hi->IonStg-1] ); 03407 03408 /* print wavelengths, either line in main printout labels, 03409 * or in various units in exponential notation - prt_wl is in prt.c */ 03410 if( strcmp( punch.chConPunEnr[punch.ipConPun], "labl" )== 0 ) 03411 { 03412 prt_wl( ioPUN , t->WLAng ); 03413 } 03414 else 03415 { 03416 /* this converts energy in Rydbergs into any of the other units */ 03417 fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t->EnergyWN * WAVNRYD)) ); 03418 } 03419 /* print the optical depth */ 03420 fprintf( ioPUN , "\t%.3f", t->Emis->TauIn ); 03421 /* damping constant */ 03422 fprintf(ioPUN, "\t%.3e", 03423 t->Emis->dampXvel / DoppVel.doppler[ t->Hi->nelem -1 ] ); 03424 fprintf(ioPUN, "\n"); 03425 } 03426 } 03427 else if( lgPopsFirstCall ) 03428 { 03429 char chLbl[11]; 03430 /* first call to line populations, print atomic parameters and indices */ 03431 strcpy( chLbl, chLineLbl(t) ); 03432 fprintf(ioPUN, "%li\t%s" , index , chLbl ); 03433 /* stat weights */ 03434 fprintf(ioPUN, "\t%.0f\t%.0f", 03435 t->Lo->g ,t->Hi->g); 03436 /* energy difference, gf */ 03437 fprintf(ioPUN, "\t%.2f\t%.3e", 03438 t->EnergyWN ,t->Emis->gf); 03439 fprintf(ioPUN, "\n"); 03440 } 03441 else 03442 { 03443 /* not first call, so do level populations and indices defined above */ 03444 if( t->Hi->Pop > xLimit ) 03445 { 03446 fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index, 03447 /* >>chng 05 may 08, add abundances, which for iso-seq species is 03448 * the density of the parent ion, for other lines, is unity. 03449 * had not been included so pops for iso seq were rel to parent ion. 03450 * caught by John Evrett */ 03451 t->Lo->Pop*abundance , 03452 t->Hi->Pop*abundance ); 03453 } 03454 } 03455 } 03456 03457 /*PunchNewContinuum produce the 'punch new continuum' output */ 03458 STATIC void PunchNewContinuum(FILE * ioPUN, long ipCon ) 03459 { 03460 long int ipLo, ipHi, 03461 j , 03462 ncells; 03463 03464 double 03465 wllo=3500. , 03466 wlhi=7000. , 03467 resolution = 2.; 03468 03469 double *energy, 03470 *cont_incid, 03471 *cont_atten, 03472 *diffuse_in, 03473 *diffuse_out, 03474 *emis_lines_out, 03475 *emis_lines_in; 03476 03477 /* get the low limit, cp_range_min is zero if not set */ 03478 wllo = MAX2( rfield.anu[0] , punch.cp_range_min[ipCon] ); 03479 if( punch.cp_range_max[ipCon] > 0. ) 03480 { 03481 /* set to smaller of these two */ 03482 wlhi = MIN2( rfield.anu[rfield.nflux-1] , punch.cp_range_max[ipCon] ); 03483 } 03484 else 03485 { 03486 /* use high-energy limit since not set */ 03487 wlhi = rfield.anu[rfield.nflux-1]; 03488 } 03489 03490 if( punch.cp_resolving_power[ipCon] != 0. ) 03491 { 03492 /* next two bogus to fool lint - they cannot happen */ 03493 ipLo = LONG_MAX; 03494 ipHi = LONG_MAX; 03495 /* we will set a new continuum mesh */ 03496 ncells = (long)((wlhi-wllo)/resolution + 1); 03497 } 03498 else 03499 { 03500 /* use native continuum mesh */ 03501 ipLo = ipoint(wllo)-1; 03502 ipHi = ipoint(wlhi)-1; 03503 ncells = ipHi - ipLo + 1; 03504 } 03505 03506 /* now allocate the space */ 03507 energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); 03508 cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); 03509 cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); 03510 diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); 03511 diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); 03512 emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); 03513 emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ); 03514 /*emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double)); 03515 emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));*/ 03516 03517 /* was the resolution changed from the default native resolution ? */ 03518 if( punch.cp_resolving_power[ipCon] != 0. ) 03519 { 03520 /* now create energy grid */ 03521 energy[0] = wlhi; 03522 j = 0; 03523 while( energy[j]-resolution > wllo ) 03524 { 03525 ++j; 03526 ASSERT( j< ncells+1 ); 03527 energy[j] = energy[j-1] - resolution; 03528 } 03529 03530 ncells = j; 03531 /* now convert to rydbergs */ 03532 for( j=0; j<ncells; ++j ) 03533 { 03534 energy[j] = RYDLAM / energy[j]; 03535 } 03536 } 03537 else 03538 { 03539 /* now convert to rydbergs */ 03540 for( j=0; j<ncells; ++j ) 03541 { 03542 energy[j] = rfield.AnuOrg[j+ipLo] - rfield.widflx[j+ipLo]/2.; 03543 } 03544 } 03545 03546 /* for cdSPEC the energy vector is the lower edge of the energy cell */ 03547 /* get incident continuum */ 03548 cdSPEC( 1 , energy , ncells , cont_incid ); 03549 /* get attenuated incident continuum */ 03550 cdSPEC( 2 , energy , ncells , cont_atten ); 03551 /* get attenuated incident continuum */ 03552 cdSPEC( 5 , energy , ncells , diffuse_in ); 03553 /* get attenuated incident continuum */ 03554 cdSPEC( 4 , energy , ncells , diffuse_out ); 03556 /* different parts of lines */ 03557 cdSPEC( 6 , energy , ncells , emis_lines_out ); 03558 cdSPEC( 7 , energy , ncells , emis_lines_in ); 03559 /*cdSPEC( 8 , energy , ncells , emis_lines_pump_out ); 03560 cdSPEC( 9 , energy , ncells , emis_lines_pump_in );*/ 03561 03562 /* for this example we will do a wavelength range */ 03563 for( j=0; j<ncells-1; ++j ) 03564 { 03565 /* photon energy in appropriate energy or wavelength units */ 03566 fprintf( ioPUN,"%.3e\t", AnuUnit((realnum)(energy[j]+rfield.widflx[j+ipLo]/2.) ) ); 03567 fprintf( ioPUN,"%.3e\t", cont_incid[j] ); 03568 fprintf( ioPUN,"%.3e\t", cont_atten[j] ); 03569 fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] ); 03570 fprintf( ioPUN,"%.3e", 03571 emis_lines_out[j]+emis_lines_in[j]/*+emis_lines_pump_out[j]+emis_lines_pump_in[j]*/ ); 03572 fprintf( ioPUN,"\n" ); 03573 } 03574 03575 free(energy); 03576 free(cont_incid); 03577 free(diffuse_in); 03578 free(diffuse_out); 03579 free(cont_atten); 03580 free(emis_lines_out); 03581 free(emis_lines_in); 03582 /*free(emis_lines_pump_out); 03583 free(emis_lines_pump_in );*/ 03584 } 03585 03586 /* punch AGN3 hemiss, for Chapter 4, routine is below */ 03587 STATIC void AGN_Hemis(FILE *ioPUN ) 03588 { 03589 const int NTE = 4; 03590 realnum te[NTE] = {5000., 10000., 15000., 20000.}; 03591 realnum *agn_continuum[NTE]; 03592 double TempSave = phycon.te; 03593 long i , j; 03594 03595 DEBUG_ENTRY( "AGN_Hemis()" ); 03596 03597 /* make table of continuous emission at various temperatuers */ 03598 /* first allocate space */ 03599 for( i=0;i<NTE; ++i) 03600 { 03601 agn_continuum[i] = (realnum *)MALLOC((unsigned)rfield.nflux*sizeof(realnum) ); 03602 03603 /* set the next temperature */ 03604 /* recompute everything at this new temp */ 03605 TempChange(te[i] , true); 03606 /* converge the pressure-temperature-ionization solution for this zone */ 03607 ConvPresTempEdenIoniz(); 03608 03609 /* now get the thermal emission */ 03610 RT_diffuse(); 03611 for(j=0;j<rfield.nflux; ++j ) 03612 { 03613 agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(realnum)dense.eden/ 03614 (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] ); 03615 } 03616 } 03617 03618 /* print title for line */ 03619 fprintf(ioPUN,"wl"); 03620 for( i=0;i<NTE; ++i) 03621 { 03622 fprintf(ioPUN,"\tT=%.0f",te[i]); 03623 } 03624 fprintf( ioPUN , "\tcont\n"); 03625 03626 /* not print all n temperatures across a line */ 03627 for(j=0;j<rfield.nflux; ++j ) 03628 { 03629 fprintf( ioPUN , "%12.5e", 03630 AnuUnit(rfield.AnuOrg[j]) ); 03631 /* loop over the temperatures, and for each, calculate a continuum */ 03632 for( i=0;i<NTE; ++i) 03633 { 03634 fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]); 03635 } 03636 /* cont label and end of line*/ 03637 fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j]); 03638 } 03639 03640 /* now free the continua */ 03641 for( i=0;i<NTE; ++i) 03642 { 03643 free( agn_continuum[i] ); 03644 } 03645 03646 /* Restore temperature stored before this routine was called */ 03647 /* and force update */ 03648 TempChange(TempSave , true); 03649 03650 fprintf( ioQQQ, "AGN_Hemis - result of punch AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n"); 03651 cdEXIT(EXIT_FAILURE); 03652 } 03653 03654 /*punResults punch results from punch results command */ 03655 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */ 03656 STATIC void punResults(FILE* ioPUN) 03657 { 03658 char chCard[INPUT_LINE_LENGTH]; 03659 bool lgEOF; 03660 long int i , nelem , ion; 03661 03662 DEBUG_ENTRY( "punResults()" ); 03663 03664 /* used to punch out line intensities, optical depths, 03665 * and column densities */ 03666 lgEOF = false; 03667 /* first initialize the line image reader */ 03668 input_init(); 03669 03670 fprintf( ioPUN, "**********************************************************************************************************************************\n" ); 03671 lgEOF = false; 03672 input_readarray(chCard,&lgEOF); 03673 while( !lgEOF ) 03674 { 03675 /* will get copy of input line image as all capital letters here */ 03676 char chCAPS[INPUT_LINE_LENGTH]; 03677 strcpy( chCAPS , chCard ); 03678 caps( chCAPS ); 03679 /* keyword HIDE means to hide the command - do not print it */ 03680 if( !nMatch( "HIDE" , chCAPS ) ) 03681 fprintf( ioPUN, "%s\n", chCard ); 03682 input_readarray(chCard,&lgEOF); 03683 } 03684 03685 /* first print any cautions or warnings */ 03686 cdWarnings(ioPUN); 03687 cdCautions(ioPUN); 03688 fprintf( ioPUN, "**********************************************************************************************************************************\n" ); 03689 03690 fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec ); 03691 03692 fprintf( ioPUN, "BEGIN EMISSION LINES\n" ); 03693 PunResults1Line(ioPUN," ",0,0.,"Start"); 03694 03695 for( i=0; i < LineSave.nsum; i++ ) 03696 { 03697 if( LineSv[i].sumlin[LineSave.lgLineEmergent] > 0. ) 03698 { 03699 PunResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength, 03700 LineSv[i].sumlin[LineSave.lgLineEmergent],"Line "); 03701 } 03702 } 03703 03704 PunResults1Line(ioPUN," ",0,0.,"Flush"); 03705 03706 fprintf( ioPUN, " \n" ); 03707 03708 fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" ); 03709 03710 /* this dumps out the whole array,*/ 03711 /* following loop relies on LIMELM being 30, assert it here in case 03712 * this is ever changed */ 03713 ASSERT( LIMELM == 30 ); 03714 /* this order of indices is to keep 30 as the fastest variable, 03715 * and the 32 (LIMELM+1) as the slower one */ 03716 for( nelem=0; nelem<LIMELM; nelem++ ) 03717 { 03718 for(ion=0; ion < nelem+1; ion++) 03719 { 03720 fprintf( ioPUN, " %10.3e", mean.xIonMeans[0][nelem][ion] ); 03721 /* throw line feed every 10 numbers */ 03722 if( nelem==9|| nelem==19 || nelem==29 ) 03723 { 03724 fprintf( ioPUN, "\n" ); 03725 } 03726 } 03727 } 03728 03729 fprintf( ioPUN, "END OF RESULTS\n" ); 03730 fprintf( ioPUN, "**********************************************************************************************************************************\n" ); 03731 return; 03732 } 03733 03734 static const int LINEWIDTH = 6; 03735 03736 /*PunResults1Line do single line of output for the punch results and punch line intensity commands */ 03737 /* the number of emission lines across one line of printout */ 03738 STATIC void PunResults1Line( 03739 FILE * ioPUN, 03740 /* 4 char + null string */ 03741 const char *chLab, 03742 realnum wl, 03743 double xInten, 03744 /* null terminated string saying what to do */ 03745 const char *chFunction) 03746 { 03747 03748 long int i; 03749 static realnum wavelength[LINEWIDTH]; 03750 static long ipLine; 03751 static double xIntenSave[LINEWIDTH]; 03752 static char chLabSave[LINEWIDTH][5]; 03753 03754 DEBUG_ENTRY( "PunResults1Line()" ); 03755 03756 /* if LineWidth is changed then change format in write too */ 03757 03758 if( strcmp(chFunction,"Start") == 0 ) 03759 { 03760 ipLine = 0; 03761 } 03762 else if( strcmp(chFunction,"Line ") == 0 ) 03763 { 03764 /* save results in array so that they can be printed when done */ 03765 wavelength[ipLine] = wl; 03766 strcpy( chLabSave[ipLine], chLab ); 03767 xIntenSave[ipLine] = xInten; 03768 03769 /* now increment the counter and then check if we have filled the line, 03770 * and so should print it */ 03771 ++ipLine; 03772 /* do print now if we are in column mode (one line per line) or if we have filled up 03773 * the line */ 03774 if( ( strcmp(punch.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH ) 03775 { 03776 /* "array " is usual array 6 wide, "column" is one line per line */ 03777 for( i=0; i < ipLine; i++ ) 03778 { 03779 fprintf( ioPUN, " %4.4s ", chLabSave[i] ); 03780 prt_wl( ioPUN, wavelength[i] ); 03781 fprintf( ioPUN,"\t%.3e", xIntenSave[i]); 03782 /* >>chng 02 apr 24, do not print type */ 03783 /* single column for input into data base */ 03784 if( strcmp(punch.chPunRltType,"column") == 0 ) 03785 fprintf( ioPUN, "\n" ); 03786 } 03787 /* only put cr if we did not just put one */ 03788 if( strcmp(punch.chPunRltType,"array ") == 0 ) 03789 fprintf( ioPUN, " \n" ); 03790 ipLine = 0; 03791 } 03792 } 03793 else if( strcmp(chFunction,"Flush") == 0 ) 03794 { 03795 if( ipLine > 0 ) 03796 { 03797 /* this is an option to print many emission lines across an output line, 03798 * the array option, or a single column of numbers, the column option 03799 * that appears on the "punch results" and "punch intensity" commands 03800 */ 03801 /* usual array 6 wide */ 03802 for( i=0; i < ipLine; i++ ) 03803 { 03804 fprintf( ioPUN, " %4.4s", chLabSave[i] ); 03805 prt_wl( ioPUN, wavelength[i] ); 03806 fprintf( ioPUN,"\t%.3e", xIntenSave[i]); 03807 /* >>chng 02 apr 24, do not print type */ 03808 /* single column for input into data base */ 03809 if( strcmp(punch.chPunRltType,"column") == 0 ) 03810 fprintf( ioPUN, "\n" ); 03811 } 03812 if( strcmp(punch.chPunRltType,"array ") == 0 ) 03813 fprintf( ioPUN, " \n" ); 03814 } 03815 } 03816 else 03817 { 03818 fprintf( ioQQQ, " PunResults1Line called with insane request=%5.5s\n", 03819 chFunction ); 03820 cdEXIT(EXIT_FAILURE); 03821 } 03822 return; 03823 } 03824 03825 static const int NENR_GAUNT = 37; 03826 static const int NTE_GAUNT = 21; 03827 03828 /*PunchGaunts called by punch gaunts command to output gaunt factors */ 03829 STATIC void PunchGaunts(FILE* ioPUN) 03830 { 03831 long int i, 03832 charge, 03833 ite, 03834 j; 03835 03836 realnum ener[NENR_GAUNT], 03837 ste[NTE_GAUNT], 03838 z; 03839 double tempsave; 03840 double g[NENR_GAUNT][NTE_GAUNT], gfac; 03841 03842 03843 DEBUG_ENTRY( "PunchGaunts()" ); 03844 03845 /* this routine is called from the PUNCH GAUNTS command 03846 * it drives the gaunt factor routine to punch gaunts over full range */ 03847 tempsave = phycon.te; 03848 03849 for( i=0; i < NTE_GAUNT; i++ ) 03850 { 03851 ste[i] = 0.5f*i; 03852 } 03853 03854 for( i=0; i < NENR_GAUNT; i++ ) 03855 { 03856 ener[i] = 0.5f*i - 8.f; 03857 } 03858 03859 for( charge = 1; charge<LIMELM; charge++ ) 03860 { 03861 /* nuclear charge */ 03862 z = (realnum)log10((double)charge); 03863 03864 /* energy is log of energy */ 03865 for( ite=0; ite < NTE_GAUNT; ite++ ) 03866 { 03867 phycon.alogte = ste[ite]; 03868 phycon.te = pow(10.,phycon.alogte); 03869 03870 for( j=0; j < NENR_GAUNT; j++ ) 03871 { 03872 gfac = cont_gaunt_calc( phycon.te, (double)charge, pow(10.,(double)ener[j]) ); 03873 /*fprintf(ioQQQ, "z %.2e ener %.2e temp %.2e gfac %.3e \n", 03874 pow(10.,z), pow(10.,ener[j]), (double)phycon.te, gfac );*/ 03875 g[j][ite] = gfac; 03876 } 03877 } 03878 03879 /* now punch out the results */ 03880 fprintf( ioPUN, " " ); 03881 for( i=1; i <= NTE_GAUNT; i++ ) 03882 { 03883 fprintf( ioPUN, "\t%6.1f", ste[i-1] ); 03884 } 03885 fprintf( ioPUN, "\n" ); 03886 03887 for( j=0; j < NENR_GAUNT; j++ ) 03888 { 03889 fprintf( ioPUN, "\t%6.2f", ener[j] ); 03890 for( ite=0; ite < NTE_GAUNT; ite++ ) 03891 { 03892 fprintf( ioPUN, "\t%6.2f", log10(g[j][ite]) ); 03893 } 03894 fprintf( ioPUN, "\n" ); 03895 } 03896 03897 fprintf( ioPUN, " " ); 03898 for( i=0; i < NTE_GAUNT; i++ ) 03899 { 03900 fprintf( ioPUN, "\t%6.1f", ste[i] ); 03901 } 03902 fprintf( ioPUN, "\n" ); 03903 03904 /* now do the same thing as triplets */ 03905 fprintf( ioPUN, " " ); 03906 for( i=0; i < NTE_GAUNT; i++ ) 03907 { 03908 fprintf( ioPUN, "\t%6.1f", ste[i] ); 03909 } 03910 fprintf( ioPUN, "\n" ); 03911 03912 for( i=0; i < NTE_GAUNT; i++ ) 03913 { 03914 for( j=0; j < NENR_GAUNT; j++ ) 03915 { 03916 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", ste[i], ener[j], 03917 log10(g[j][i]) ); 03918 } 03919 } 03920 03921 fprintf( ioPUN, "Below is log(u), log(gamma**2), gff\n" ); 03922 /* do the same thing as above, but this time print log(u) and log(gamma2) instead of temp and energy. */ 03923 for( i=0; i < NTE_GAUNT; i++ ) 03924 { 03925 for( j=0; j < NENR_GAUNT; j++ ) 03926 { 03927 fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", 2.*z + log10(TE1RYD) - ste[i] , log10(TE1RYD)+ener[j]-ste[i], 03928 log10(g[j][i]) ); 03929 } 03930 } 03931 fprintf( ioPUN, "end of charge = %li\n", charge); 03932 fprintf( ioPUN, "****************************\n"); 03933 } 03934 03935 phycon.te = tempsave; 03936 phycon.alogte = log10( phycon.te ); 03937 return; 03938 }