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 /*PrtZone print out individual zone results */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "iso.h" 00007 #include "grainvar.h" 00008 #include "pressure.h" 00009 #include "wind.h" 00010 #include "conv.h" 00011 #include "trace.h" 00012 #include "magnetic.h" 00013 #include "dense.h" 00014 #include "called.h" 00015 #include "dynamics.h" 00016 #include "h2.h" 00017 #include "mole.h" 00018 #include "secondaries.h" 00019 #include "opacity.h" 00020 #include "colden.h" 00021 #include "geometry.h" 00022 #include "hmi.h" 00023 #include "rfield.h" 00024 #include "thermal.h" 00025 #include "radius.h" 00026 #include "phycon.h" 00027 #include "abund.h" 00028 #include "hydrogenic.h" 00029 #include "ionbal.h" 00030 #include "elementnames.h" 00031 #include "atomfeii.h" 00032 #include "prt.h" 00033 00034 void PrtZone(void) 00035 { 00036 char chField7[8]; 00037 char chLet, 00038 chQHMark; 00039 long int i, 00040 ishift, 00041 nelem , 00042 nd, 00043 mol; 00044 double cdif, 00045 coninc, 00046 con_density, 00047 fac, 00048 hatmic; 00049 00050 DEBUG_ENTRY( "PrtZone()" ); 00051 00052 if( thermal.lgUnstable ) 00053 { 00054 chLet = 'u'; 00055 } 00056 else 00057 { 00058 chLet = ' '; 00059 } 00060 00061 /* middle of zone for printing 00062 rmidle = radius.Radius - radius.drad*0.5*radius.dRadSign; 00063 dmidle = radius.depth - radius.drad*0.5; */ 00064 00065 /* option to print single line when quiet but tracing convergence 00066 * with "trace convergence" command */ 00067 if( called.lgTalk || trace.nTrConvg ) 00068 { 00069 /* print either ####123 or ###1234 */ 00070 if( nzone <= 999 ) 00071 { 00072 sprintf( chField7, "####%3ld", nzone ); 00073 } 00074 else 00075 { 00076 sprintf( chField7, "###%4ld", nzone ); 00077 } 00078 00079 fprintf(ioQQQ, " %7.7s %cTe:",chField7, chLet); 00080 PrintE93(ioQQQ,phycon.te); 00081 fprintf(ioQQQ," Hden:"); 00082 PrintE93(ioQQQ,dense.gas_phase[ipHYDROGEN]); 00083 fprintf(ioQQQ," Ne:"); 00084 PrintE93(ioQQQ,dense.eden); 00085 fprintf(ioQQQ," R:"); 00086 PrintE93(ioQQQ,radius.Radius_mid_zone ); 00087 fprintf(ioQQQ," R-R0:"); 00088 PrintE93(ioQQQ,radius.depth_mid_zone); 00089 fprintf(ioQQQ," dR:"); 00090 PrintE93(ioQQQ,radius.drad); 00091 fprintf(ioQQQ," NTR:%3ld Htot:",conv.nPres2Ioniz); 00092 PrintE93(ioQQQ,thermal.htot); 00093 fprintf(ioQQQ," T912:"); 00094 fprintf(ioQQQ,PrintEfmt("%9.2e",opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] )); 00095 fprintf(ioQQQ,"###\n"); 00096 00097 if( trace.nTrConvg ) 00098 { 00099 fprintf( ioQQQ, " H:%.2e %.2e 2H2/H: %.2e He: %.2e %.2e %.2e\n", 00100 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN], 00101 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN], 00102 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN], 00103 dense.xIonDense[ipHELIUM][0]/SDIV(dense.gas_phase[ipHELIUM]), 00104 dense.xIonDense[ipHELIUM][1]/SDIV(dense.gas_phase[ipHELIUM]), 00105 dense.xIonDense[ipHELIUM][2]/SDIV(dense.gas_phase[ipHELIUM]) 00106 ); 00107 } 00108 } 00109 00110 /* now return if not talking */ 00111 if( !called.lgTalk || trace.nTrConvg ) 00112 { 00113 return; 00114 } 00115 00116 /* lgDenFlucOn set to true in zero, only false when variable abundances are on, 00117 * lgAbTaON set true when element table used */ 00118 if( !dense.lgDenFlucOn || abund.lgAbTaON ) 00119 { 00120 fprintf( ioQQQ, " Abun:" ); 00121 for( i=0; i < LIMELM; i++ ) 00122 { 00123 fprintf( ioQQQ,PrintEfmt("%8.1e", dense.gas_phase[i] )); 00124 } 00125 fprintf( ioQQQ, "\n" ); 00126 } 00127 00128 /*------------------------------------------------- 00129 * print wind parameters if windy model */ 00130 fac = wind.windv*dense.gas_phase[ipHYDROGEN]*radius.r1r0sq; 00131 if( wind.windv != 0. ) 00132 { 00133 /* find denominator for fractional contributions */ 00134 if( wind.AccelTot == 0. ) 00135 fac = 1.; 00136 else 00137 fac = wind.AccelTot; 00138 fprintf( ioQQQ, 00139 " Dynamics wind V:%.3e km/s a(grav):%.2e a(tot):%.2e Fr(cont):%6.3f " 00140 "Fr(line):%6.3f Fr(dP):%6.3f\n", 00141 wind.windv/1e5 , 00142 -wind.AccelGravity, 00143 wind.AccelTot, 00144 wind.AccelCont/ fac, 00145 wind.AccelLine/fac, 00146 wind.AccelPres/fac ); 00147 00148 /* print advection information */ 00149 if( dynamics.lgAdvection ) 00150 DynaPrtZone(); 00151 } 00152 00153 /* print line with radiation pressure if significant */ 00154 if( pressure.pbeta > .05 ) 00155 PrtLinePres(); 00156 00157 /*---------------------------------------------------- */ 00158 00159 hatmic = 0.; 00160 for(mol = 0; mol < N_H_MOLEC; mol++) { 00161 hatmic += hmi.Hmolec[mol]*hmi.nProton[mol]; 00162 } 00163 ASSERT(hatmic > 0.); 00164 hatmic = (dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1])/hatmic; 00165 00166 fprintf( ioQQQ, " Hydrogen "); 00167 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])); 00168 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])); 00169 fprintf( ioQQQ, " H+o/Hden"); 00170 fprintf(ioQQQ,PrintEfmt("%9.2e",hatmic )); 00171 fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.Hmolec[ipMHm]/dense.gas_phase[ipHYDROGEN] )); 00172 fprintf( ioQQQ, " H- H2"); 00173 /* this is total H2, the sum of "ground" and excited */ 00174 fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.H2_total/dense.gas_phase[ipHYDROGEN])); 00175 fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.Hmolec[ipMH2p]/dense.gas_phase[ipHYDROGEN])); 00176 fprintf( ioQQQ, " H2+ HeH+"); 00177 fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.Hmolec[ipMHeHp]/dense.gas_phase[ipHYDROGEN])); 00178 fprintf( ioQQQ, " Ho+ ColD"); 00179 fprintf(ioQQQ,PrintEfmt("%9.2e",colden.colden[ipCOL_H0])); 00180 fprintf(ioQQQ,PrintEfmt("%9.2e",colden.colden[ipCOL_Hp])); 00181 fprintf( ioQQQ, "\n"); 00182 00183 /* print departure coef if desired */ 00184 if( iso.lgPrtDepartCoef[ipH_LIKE][ipHYDROGEN] ) 00185 { 00186 fprintf( ioQQQ, " Hydrogen " ); 00187 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s])); 00188 fprintf(ioQQQ,PrintEfmt("%9.2e", 1.)); 00189 fprintf( ioQQQ, " H+o/Hden"); 00190 fprintf(ioQQQ,PrintEfmt("%9.2e", (dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1])/dense.gas_phase[ipHYDROGEN])); 00191 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.hmidep)); 00192 fprintf( ioQQQ, " H- H2"); 00193 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.h2dep)); 00194 fprintf( ioQQQ, " H2+"); 00195 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.h2pdep)); 00196 fprintf( ioQQQ, " H3+"); 00197 fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.h3pdep)); 00198 fprintf( ioQQQ, "\n" ); 00199 } 00200 00201 if( prt.lgPrintHeating ) 00202 { 00203 fprintf( ioQQQ, " "); 00204 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][0]/thermal.htot)); 00205 fprintf( ioQQQ," "); 00206 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][15]/thermal.htot)); 00207 fprintf( ioQQQ," "); 00208 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][16]/thermal.htot)); 00209 fprintf( ioQQQ,"\n"); 00210 } 00211 00212 /* temperature correspoding to radiation fields */ 00213 coninc = 0.; 00214 cdif = 0.; 00215 for( i=0; i < rfield.nflux; i++ ) 00216 { 00217 /* integrated energy flux, ergs s^-1 cm^-2 */ 00218 coninc += rfield.flux[i]*(rfield.anu[i]*EN1RYD); 00219 cdif += (rfield.outlin[i] + rfield.outlin_noplot[i] + 00220 rfield.ConInterOut[i])* (rfield.anu[i]*EN1RYD); 00221 } 00222 /* convert flux in attenuated incident continuum, diffuse emission, into equivalent temperature */ 00223 coninc = pow(coninc/SPEEDLIGHT/7.56464e-15,0.25); 00224 cdif = pow(cdif/SPEEDLIGHT/7.56464e-15,0.25); 00225 00226 /* convert sum of flux into energy density, then equivalent pressure */ 00227 con_density = (coninc + cdif) / SPEEDLIGHT; 00228 con_density /= BOLTZMANN; 00229 00230 if( prt.lgPrintHeating ) 00231 { 00232 fprintf( ioQQQ, " "); 00233 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][1]/thermal.htot )); 00234 fprintf( ioQQQ, " "); 00235 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00236 fprintf( ioQQQ, " BoundCom"); 00237 fprintf(ioQQQ,PrintEfmt("%9.2e", ionbal.CompRecoilHeatLocal/ thermal.htot)); 00238 fprintf( ioQQQ, " Extra:"); 00239 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][20]/thermal.htot)); 00240 fprintf( ioQQQ, " Pairs:"); 00241 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][21]/ thermal.htot )); 00242 fprintf( ioQQQ," H-lines\n"); 00243 } 00244 00245 /* Helium */ 00246 if( dense.lgElmtOn[ipHELIUM] ) 00247 { 00248 fprintf( ioQQQ, " Helium " ); 00249 for( i=0; i < 3; i++ ) 00250 { 00251 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipHELIUM][i]/dense.gas_phase[ipHELIUM]) ); 00252 } 00253 00254 fprintf( ioQQQ, " He I2SP3"); 00255 fprintf(ioQQQ,PrintEfmt("%9.2e", 00256 StatesElem[ipHE_LIKE][ipHELIUM][ipHe2s3S].Pop*dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM] )); 00257 fprintf( ioQQQ, " Comp H,C"); 00258 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmheat )); 00259 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmcool*phycon.te)); 00260 fprintf( ioQQQ , " Fill Fac"); 00261 fprintf(ioQQQ,PrintEfmt("%9.2e", geometry.FillFac)); 00262 fprintf( ioQQQ , " Gam1/tot"); 00263 fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.H_ion_frac_photo)); 00264 fprintf( ioQQQ, "\n"); 00265 00266 /* option to print departure coef */ 00267 if( iso.lgPrtDepartCoef[ipH_LIKE][ipHELIUM] ) 00268 { 00269 fprintf( ioQQQ, " Helium " ); 00270 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipHE_LIKE][ipHELIUM][0])); 00271 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipH_LIKE][ipHELIUM][ipH1s])); 00272 fprintf(ioQQQ,PrintEfmt("%9.2e", 1.)); 00273 00274 fprintf( ioQQQ, " Comp H,C"); 00275 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmheat )); 00276 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmcool*phycon.te )); 00277 fprintf( ioQQQ , " Fill Fac"); 00278 fprintf(ioQQQ,PrintEfmt("%9.2e", geometry.FillFac )); 00279 fprintf( ioQQQ , " Gam1/tot"); 00280 fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.H_ion_frac_photo)); 00281 fprintf( ioQQQ, "\n"); 00282 } 00283 00284 /* print heating from He (and others) if desired 00285 * entry "lines" is induced line heating 00286 * 1,12 ffheat: 2,3 he triplets, 1,20 compton */ 00287 if( prt.lgPrintHeating ) 00288 { 00289 /*fprintf( ioQQQ, " %10.3e%10.3e Lines:%10.2e%10.2e Compton:%10.3e FF Heatig%10.3e\n", 00290 thermal.heating[1][0]/thermal.htot, thermal.heating[1][1]/ 00291 thermal.htot, thermal.heating[0][22]/thermal.htot, thermal.heating[1][2]/ 00292 thermal.htot, thermal.heating[0][19]/thermal.htot, thermal.heating[0][11]/ 00293 thermal.htot );*/ 00294 fprintf( ioQQQ, " "); 00295 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][0]/thermal.htot)); 00296 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][1]/thermal.htot)); 00297 fprintf( ioQQQ, " Lines:"); 00298 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][22]/thermal.htot)); 00299 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][2]/thermal.htot)); 00300 fprintf( ioQQQ, " Compton:"); 00301 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][19]/thermal.htot)); 00302 fprintf( ioQQQ, " FFHeatig"); 00303 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][11]/thermal.htot)); 00304 fprintf( ioQQQ, "\n"); 00305 } 00306 00307 if( dense.lgElmtOn[ipHELIUM] ) 00308 { 00309 /* helium singlets and triplets relative to total helium gas phase density */ 00310 fac = dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM]; 00311 fprintf( ioQQQ, " He singlet n " ); 00312 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][1][ipHe1s1S].Pop*fac )); 00313 /* singlet n=2 complex */ 00314 if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]>= ipHe2p1P ) 00315 { 00316 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][1][ipHe2s1S].Pop*fac )); 00317 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][1][ipHe2p1P].Pop*fac )); 00318 } 00319 else 00320 { 00321 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00322 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00323 } 00324 /* singlet n=3 complex */ 00325 if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]>= ipHe3p1P ) 00326 { 00327 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][ipHELIUM][ipHe3s1S].Pop*fac )); 00328 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][ipHELIUM][ipHe3p1P].Pop*fac )); 00329 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][ipHELIUM][ipHe3d1D].Pop*fac )); 00330 } 00331 else 00332 { 00333 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00334 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00335 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00336 } 00337 00338 fprintf( ioQQQ, " He tripl" ); 00339 /* triplet n=2 complex */ 00340 if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]>= ipHe2p3P2 ) 00341 { 00342 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][ipHELIUM][ipHe2s3S].Pop*fac )); 00343 fprintf(ioQQQ,PrintEfmt("%9.2e", 00344 StatesElem[ipHE_LIKE][ipHELIUM][ipHe2p3P0].Pop*fac+ 00345 StatesElem[ipHE_LIKE][ipHELIUM][ipHe2p3P1].Pop*fac+ 00346 StatesElem[ipHE_LIKE][ipHELIUM][ipHe2p3P2].Pop*fac )); 00347 } 00348 else 00349 { 00350 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00351 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00352 } 00353 /* triplet n=3 complex */ 00354 if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]> ipHe3d3D ) 00355 { 00356 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][ipHELIUM][ipHe3s3S].Pop*fac )); 00357 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][ipHELIUM][ipHe3p3P].Pop*fac )); 00358 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipHE_LIKE][ipHELIUM][ipHe3d3D].Pop*fac )); 00359 } 00360 else 00361 { 00362 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00363 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00364 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. )); 00365 } 00366 fprintf( ioQQQ, "\n" ); 00367 } 00368 } 00369 00370 /* loop over iso sequences to see if any populations 00371 * and/or departure coefficients need to be printed */ 00372 for( long ipISO = ipH_LIKE; ipISO < NISO; ipISO++ ) 00373 { 00374 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00375 { 00376 if( dense.lgElmtOn[nelem] ) 00377 { 00378 if( iso.lgPrtLevelPops[ipISO][nelem] ) 00379 { 00380 iso_prt_pops(ipISO, nelem, false); 00381 } 00382 if( iso.lgPrtDepartCoef[ipISO][nelem] ) 00383 { 00384 /* true says print departure coefficients 00385 * instead of populations. */ 00386 iso_prt_pops(ipISO, nelem, true); 00387 } 00388 } 00389 } 00390 } 00391 00392 /* >>chng 01 dec 08, move pressure to line before grains, after radiation properties */ 00393 /* gas pressure, pressure due to incident radiation field, rad accel */ 00394 fprintf( ioQQQ, " Pressure NgasTgas"); 00395 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresGasCurr/BOLTZMANN)); 00396 fprintf( ioQQQ, " P(total)"); 00397 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresTotlCurr)); 00398 fprintf( ioQQQ, " P( gas )"); 00399 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresGasCurr)); 00400 fprintf( ioQQQ, " P(Radtn)"); 00401 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.pres_radiation_lines_curr)); 00402 fprintf( ioQQQ, " Rad accl"); 00403 fprintf(ioQQQ,PrintEfmt("%9.2e", wind.AccelTot)); 00404 fprintf( ioQQQ, " ForceMul"); 00405 fprintf(ioQQQ,PrintEfmt("%9.2e", wind.fmul)); 00406 fprintf( ioQQQ, "\n" ); 00407 00408 fprintf( ioQQQ , " Texc(La)"); 00409 fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.TexcLya )); 00410 fprintf( ioQQQ , " T(contn)"); 00411 fprintf(ioQQQ,PrintEfmt("%9.2e", coninc )); 00412 fprintf( ioQQQ , " T(diffs)"); 00413 fprintf(ioQQQ,PrintEfmt("%9.2e", cdif )); 00414 /* print the total radiation density expressed as an equivalent gas pressure */ 00415 fprintf( ioQQQ , " nT (c+d)"); 00416 fprintf(ioQQQ,PrintEfmt("%9.2e", con_density )); 00417 /* print the radiation to gas pressure */ 00418 fprintf( ioQQQ , " Prad/Gas"); 00419 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.pbeta )); 00420 /* magnetic to gas pressure ratio */ 00421 fprintf( ioQQQ , " Pmag/Gas"); 00422 fprintf(ioQQQ,PrintEfmt("%9.2e", magnetic.pressure / pressure.PresGasCurr) ); 00423 fprintf( ioQQQ, "\n" ); 00424 00425 if( gv.lgGrainPhysicsOn ) 00426 { 00427 for( nd=0; nd < gv.nBin; nd++ ) 00428 { 00429 /* Change things so the quantum heated dust species are marked with an 00430 * asterisk just after the name (K Volk) 00431 * added QHMARK here and in the write statement */ 00432 chQHMark = (char)(( gv.bin[nd]->lgQHeat && gv.bin[nd]->lgUseQHeat ) ? '*' : ' '); 00433 fprintf( ioQQQ, "%-12.12s%c DustTemp",gv.bin[nd]->chDstLab, chQHMark); 00434 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->tedust)); 00435 fprintf( ioQQQ, " Pot Volt"); 00436 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->dstpot*EVRYD)); 00437 fprintf( ioQQQ, " Chrg (e)"); 00438 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->AveDustZ)); 00439 fprintf( ioQQQ, " drf cm/s"); 00440 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->DustDftVel)); 00441 fprintf( ioQQQ, " Heating:"); 00442 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->GasHeatPhotoEl)); 00443 fprintf( ioQQQ, " Frac tot"); 00444 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->GasHeatPhotoEl/thermal.htot)); 00445 fprintf( ioQQQ, "\n" ); 00446 } 00447 } 00448 /* >>chng 00 apr 20, moved punch-out of quantum heating data to qheat(), by PvH */ 00449 00450 /* heavy element molecules */ 00451 if( findspecies("CO")->hevmol > 0. ) 00452 { 00453 fprintf( ioQQQ, " Molecules CH/Ctot:"); 00454 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("CH")->hevmol/dense.gas_phase[ipCARBON])); 00455 fprintf( ioQQQ, " CH+/Ctot"); 00456 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("CH+")->hevmol/dense.gas_phase[ipCARBON])); 00457 fprintf( ioQQQ, " CO/Ctot:"); 00458 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("CO")->hevmol/dense.gas_phase[ipCARBON])); 00459 fprintf( ioQQQ, " CO+/Ctot"); 00460 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("CO+")->hevmol/dense.gas_phase[ipCARBON])); 00461 fprintf( ioQQQ, " H2O/Otot"); 00462 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("H2O")->hevmol/dense.gas_phase[ipOXYGEN])); 00463 fprintf( ioQQQ, " OH/Ototl"); 00464 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("OH")->hevmol/dense.gas_phase[ipOXYGEN])); 00465 fprintf( ioQQQ, "\n"); 00466 } 00467 00468 /* information about the large H2 molecule - this just returns if not turned on */ 00469 H2_Prt_Zone(); 00470 00471 /* Lithium, Beryllium */ 00472 if( dense.lgElmtOn[ipLITHIUM] || dense.lgElmtOn[ipBERYLLIUM] || 00473 (secondaries.csupra[ipHYDROGEN][0]>0.) ) 00474 { 00475 fprintf( ioQQQ, " Lithium " ); 00476 for( i=0; i < 4; i++ ) 00477 { 00478 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipLITHIUM][i]/MAX2(1e-35,dense.gas_phase[ipLITHIUM]) )); 00479 } 00480 fprintf( ioQQQ, " Berylliu" ); 00481 for( i=0; i < 5; i++ ) 00482 { 00483 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipBERYLLIUM][i]/MAX2(1e-35,dense.gas_phase[ipBERYLLIUM])) ); 00484 } 00485 00486 /* print secondary ionization rate for atomic hydrogen */ 00487 fprintf( ioQQQ, " sec ion:" ); 00488 fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.csupra[ipHYDROGEN][0]) ); 00489 fprintf( ioQQQ, "\n" ); 00490 00491 /* option to print heating due to these stages*/ 00492 if( prt.lgPrintHeating ) 00493 { 00494 fprintf( ioQQQ, " " ); 00495 for( i=0; i < 3; i++ ) 00496 { 00497 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipLITHIUM][i]/ thermal.htot) ); 00498 } 00499 fprintf( ioQQQ, " " ); 00500 00501 for( i=0; i < 4; i++ ) 00502 { 00503 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipBERYLLIUM][i]/thermal.htot )); 00504 } 00505 fprintf( ioQQQ, "\n" ); 00506 } 00507 } 00508 00509 /* Boron */ 00510 if( dense.lgElmtOn[ipBORON] ) 00511 { 00512 fprintf( ioQQQ, " Boron " ); 00513 for( i=0; i < 6; i++ ) 00514 { 00515 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipBORON][i]/MAX2(1e-35,dense.gas_phase[ipBORON]) )); 00516 } 00517 fprintf( ioQQQ, "\n" ); 00518 00519 /* option to print heating*/ 00520 if( prt.lgPrintHeating ) 00521 { 00522 fprintf( ioQQQ, " " ); 00523 for( i=0; i < 5; i++ ) 00524 { 00525 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipBORON][i]/thermal.htot )); 00526 } 00527 fprintf( ioQQQ, "\n" ); 00528 } 00529 } 00530 00531 /* Carbon */ 00532 fprintf( ioQQQ, " Carbon " ); 00533 for( i=0; i < 7; i++ ) 00534 { 00535 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipCARBON][i]/SDIV(dense.gas_phase[ipCARBON])) ); 00536 } 00537 /* some molecules trail the line */ 00538 fprintf( ioQQQ, " H2O+/O " ); 00539 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("H2O+")->hevmol/MAX2(1e-35,dense.gas_phase[ipOXYGEN]) )); 00540 fprintf( ioQQQ, " OH+/Otot" ); 00541 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("OH+")->hevmol/ MAX2(1e-35,dense.gas_phase[ipOXYGEN]) )); 00542 /* print extra heating, normally zero */ 00543 fprintf( ioQQQ, " Hex(tot)" ); 00544 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][20] )); 00545 fprintf( ioQQQ, "\n" ); 00546 00547 /* option to print heating*/ 00548 if( prt.lgPrintHeating ) 00549 { 00550 fprintf( ioQQQ, " " ); 00551 for( i=0; i < ipCARBON+1; i++ ) 00552 { 00553 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipCARBON][i]/ thermal.htot) ); 00554 } 00555 fprintf( ioQQQ, "\n" ); 00556 } 00557 00558 /* Nitrogen */ 00559 fprintf( ioQQQ, " Nitrogen " ); 00560 for( i=1; i <= 8; i++ ) 00561 { 00562 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipNITROGEN][i-1]/ SDIV(dense.gas_phase[ipNITROGEN]) )); 00563 } 00564 fprintf( ioQQQ, " O2/Ototl" ); 00565 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("O2")->hevmol/MAX2(1e-35,dense.gas_phase[ipOXYGEN]))); 00566 fprintf( ioQQQ, " O2+/Otot" ); 00567 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecies("O2+")->hevmol/ MAX2(1e-35,dense.gas_phase[ipOXYGEN]) )); 00568 fprintf( ioQQQ, "\n" ); 00569 00570 /* option to print heating*/ 00571 if( prt.lgPrintHeating ) 00572 { 00573 fprintf( ioQQQ, " " ); 00574 for( i=0; i < ipNITROGEN+1; i++ ) 00575 { 00576 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipNITROGEN][i]/ thermal.htot )); 00577 } 00578 fprintf( ioQQQ, "\n" ); 00579 } 00580 00581 # if 0 00582 /* Oxygen */ 00583 fprintf( ioQQQ, " Oxygen " ); 00584 for( i=1; i <= 9; i++ ) 00585 { 00586 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipOXYGEN][i-1]/ SDIV(dense.gas_phase[ipOXYGEN]) )); 00587 } 00588 fprintf( ioQQQ, "\n" ); 00589 00590 /* option to print heating*/ 00591 if( prt.lgPrintHeating ) 00592 { 00593 fprintf( ioQQQ, " " ); 00594 for( i=0; i < ipOXYGEN+1; i++ ) 00595 { 00596 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipOXYGEN][i]/ thermal.htot )); 00597 } 00598 fprintf( ioQQQ, "\n" ); 00599 } 00600 # endif 00601 00602 /* now print rest of elements inside loops */ 00603 /* fluorine through Magnesium */ 00604 for( nelem=ipOXYGEN; nelem < ipALUMINIUM; ++nelem ) 00605 { 00606 if( dense.lgElmtOn[nelem] ) 00607 { 00608 /* print the element name and amount of shift */ 00609 fprintf( ioQQQ, " %10.10s ", elementnames.chElementName[nelem]); 00610 00611 for( i=0; i < nelem+2; i++ ) 00612 { 00613 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[nelem][i]/dense.gas_phase[nelem] )); 00614 } 00615 fprintf( ioQQQ, "\n" ); 00616 00617 /* print heating but only if needed */ 00618 if( prt.lgPrintHeating ) 00619 { 00620 fprintf( ioQQQ, " " ); 00621 for( i=0; i < nelem+1; i++ ) 00622 { 00623 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[nelem][i]/thermal.htot )); 00624 } 00625 fprintf( ioQQQ, "\n" ); 00626 } 00627 } 00628 } 00629 00630 /* Aluminium through Zinc */ 00631 for( nelem=ipALUMINIUM; nelem < LIMELM; ++nelem ) 00632 { 00633 if( dense.lgElmtOn[nelem] ) 00634 { 00635 /* number of ionization stages to print across the page */ 00636 /*@-redef@*/ 00637 enum {LINE=13}; 00638 /*@+redef@*/ 00639 ishift = MAX2(0,dense.IonHigh[nelem]-LINE+1); 00640 00641 /* print the element name and amount of shift */ 00642 fprintf( ioQQQ, " %10.10s%2ld ", elementnames.chElementName[nelem],ishift ); 00643 00644 for( i=0; i < LINE; i++ ) 00645 { 00646 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[nelem][i+ishift]/dense.gas_phase[nelem]) ); 00647 } 00648 fprintf( ioQQQ, "\n" ); 00649 00650 /* print heating but only if needed */ 00651 if( prt.lgPrintHeating ) 00652 { 00653 fprintf( ioQQQ, " " ); 00654 for( i=0; i < LINE; i++ ) 00655 { 00656 fprintf(ioQQQ, 00657 PrintEfmt("%9.2e", thermal.heating[nelem][i+ishift]/thermal.htot )); 00658 } 00659 fprintf( ioQQQ, "\n" ); 00660 } 00661 } 00662 } 00663 00664 /* call FeII print routine if large atom is turned on */ 00665 FeIIPrint(); 00666 return; 00667 } 00668 00669