00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "iso.h"
00009 #include "cddrive.h"
00010 #include "physconst.h"
00011 #include "lines.h"
00012 #include "taulines.h"
00013 #include "warnings.h"
00014 #include "phycon.h"
00015 #include "dense.h"
00016 #include "grainvar.h"
00017 #include "h2.h"
00018 #include "hmi.h"
00019 #include "thermal.h"
00020 #include "hydrogenic.h"
00021 #include "rt.h"
00022 #include "atmdat.h"
00023 #include "timesc.h"
00024 #include "opacity.h"
00025 #include "struc.h"
00026 #include "pressure.h"
00027 #include "conv.h"
00028 #include "geometry.h"
00029 #include "called.h"
00030 #include "iterations.h"
00031 #include "version.h"
00032 #include "colden.h"
00033 #include "input.h"
00034 #include "rfield.h"
00035 #include "radius.h"
00036 #include "peimbt.h"
00037 #include "oxy.h"
00038 #include "ipoint.h"
00039 #include "lines_service.h"
00040 #include "mean.h"
00041 #include "wind.h"
00042 #include "prt.h"
00043
00044
00045 STATIC void gett2o3(realnum *tsqr);
00046
00047
00048 STATIC void gett2(realnum *tsqr);
00049
00050
00051 void PrtFinal(void)
00052 {
00053 short int *lgPrt;
00054 realnum *wavelength;
00055 realnum *sclsav , *scaled;
00056 long int *ipSortLines;
00057 double *xLog_line_lumin;
00058 char **chSLab;
00059 char *chSTyp;
00060 char chCKey[5],
00061 chGeo[7],
00062 chPlaw[21];
00063
00064 long int
00065 i,
00066 ipEmType ,
00067 ipNegIntensity[33],
00068 ip2500,
00069 ip2kev,
00070 iprnt,
00071 j,
00072 nd,
00073 nline,
00074 nNegIntenLines;
00075 double o4363,
00076 bacobs,
00077 absint,
00078 bacthn,
00079 hbcab,
00080 hbeta,
00081 o5007;
00082
00083 double a,
00084 ajmass,
00085 ajmin,
00086 alfox,
00087 bot,
00088 bottom,
00089 bremtk,
00090 coleff,
00091
00092 d[8],
00093 dmean,
00094 ferode,
00095 he4686,
00096 he5876,
00097 heabun,
00098 hescal,
00099 pion,
00100 plow,
00101 powerl,
00102 qion,
00103 qlow,
00104 rate,
00105 ratio,
00106 reclin,
00107 rjeans,
00108 snorm,
00109 tmean,
00110 top,
00111 THI,
00112 uhel,
00113 uhl,
00114 usp,
00115 wmean,
00116 znit;
00117
00118 double bac,
00119 tel,
00120 x;
00121
00122 DEBUG_ENTRY( "PrtFinal()" );
00123
00124
00125
00126 if( !called.lgTalk || (lgAbort && nzone< 1) )
00127 {
00128 return;
00129 }
00130
00131
00132
00133
00134 ASSERT( LineSave.nsum > 1 );
00135
00136
00137 fprintf( ioQQQ, "\f\n");
00138 fprintf( ioQQQ, "%23c", ' ' );
00139 int len = (int)strlen(t_version::Inst().chVersion);
00140 int repeat = (72-len)/2;
00141 for( i=0; i < repeat; ++i )
00142 fprintf( ioQQQ, "*" );
00143 fprintf( ioQQQ, "> Cloudy %s <", t_version::Inst().chVersion );
00144 for( i=0; i < 72-repeat-len; ++i )
00145 fprintf( ioQQQ, "*" );
00146 fprintf( ioQQQ, "\n" );
00147
00148 for( i=0; i <= input.nSave; i++ )
00149 {
00150
00151
00152
00153 cap4(chCKey,input.chCardSav[i]);
00154
00155
00156 strcpy( input.chCARDCAPS , input.chCardSav[i] );
00157 caps( input.chCARDCAPS );
00158
00159
00160
00161 if( (strcmp(chCKey,"CONT")!= 0) && !nMatch( "HIDE" , input.chCARDCAPS) )
00162 fprintf( ioQQQ, "%23c* %-80s*\n", ' ', input.chCardSav[i] );
00163 }
00164
00165
00166 if( rfield.uh > 0. )
00167 {
00168 a = log10(rfield.uh);
00169 }
00170 else
00171 {
00172 a = -37.;
00173 }
00174
00175 fprintf( ioQQQ,
00176 " *********************************> Log(U):%6.2f <*********************************\n",
00177 a );
00178
00179 if( t_version::Inst().nBetaVer > 0 )
00180 {
00181 fprintf( ioQQQ,
00182 "\n This is a beta test version of the code, and is intended for testing only.\n\n" );
00183 }
00184
00185 if( warnings.lgWarngs )
00186 {
00187 fprintf( ioQQQ, " \n" );
00188 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00189 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00190 fprintf( ioQQQ, " >>>>>>>>>> Warning! Warnings exist, this calculation has serious problems.\n" );
00191 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00192 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00193 fprintf( ioQQQ, " \n" );
00194 }
00195 else if( warnings.lgCautns )
00196 {
00197 fprintf( ioQQQ,
00198 " >>>>>>>>>> Cautions are present.\n" );
00199 }
00200
00201 if( conv.lgBadStop )
00202 {
00203 fprintf( ioQQQ,
00204 " >>>>>>>>>> The calculation stopped for unintended reasons.\n" );
00205 }
00206
00207 if( iterations.lgIterAgain )
00208 {
00209 fprintf( ioQQQ,
00210 " >>>>>>>>>> Another iteration is needed.\n" );
00211 }
00212
00213
00214 if( geometry.lgSphere )
00215 {
00216 strcpy( chGeo, "Closed" );
00217 }
00218 else
00219 {
00220 strcpy( chGeo, " Open" );
00221 }
00222
00223
00224 if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
00225 {
00226 strcpy( chPlaw, " Constant Pressure " );
00227 }
00228
00229 else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
00230 {
00231 strcpy( chPlaw, " Constant Density " );
00232 }
00233
00234 else if( (strcmp(dense.chDenseLaw,"POWD") == 0 || strcmp(dense.chDenseLaw
00235 ,"POWR") == 0) || strcmp(dense.chDenseLaw,"POWC") == 0 )
00236 {
00237 strcpy( chPlaw, " Power Law Density " );
00238 }
00239
00240 else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
00241 {
00242 strcpy( chPlaw, " Rapid Fluctuations " );
00243 }
00244
00245 else if( strncmp(dense.chDenseLaw , "DLW" , 3) == 0 )
00246 {
00247 strcpy( chPlaw, " Special Density Lw " );
00248 }
00249
00250 else if( strcmp(dense.chDenseLaw,"HYDR") == 0 )
00251 {
00252 strcpy( chPlaw, " Hydrostatic Equlib " );
00253 }
00254
00255 else if( strcmp(dense.chDenseLaw,"WIND") == 0 )
00256 {
00257 strcpy( chPlaw, " Radia Driven Wind " );
00258 }
00259
00260 else if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00261 {
00262 strcpy( chPlaw, " Globule " );
00263 }
00264
00265 else
00266 {
00267 strcpy( chPlaw, " UNRECOGNIZED CPRES " );
00268 }
00269
00270 fprintf( ioQQQ,
00271 "\n Emission Line Spectrum. %20.20sModel. %6.6s geometry. Iteration %ld of %ld.\n",
00272 chPlaw, chGeo, iteration, iterations.itermx + 1 );
00273
00274
00275
00276 if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth )
00277 {
00278 fprintf( ioQQQ,
00279 " Flux observed at the Earth (erg/s/cm^2)\n" );
00280 }
00281
00282 else if( prt.lgSurfaceBrightness )
00283 {
00284 fprintf( ioQQQ,
00285 " Surface brightness (erg/s/cm^2/" );
00286 if( prt.lgSurfaceBrightness_SR )
00287 {
00288 fprintf( ioQQQ, "sr)\n");
00289 }
00290 else
00291 {
00292 fprintf( ioQQQ, "srcsec^2)\n");
00293 }
00294 }
00295
00296 else if( radius.lgPredLumin )
00297 {
00298
00299 if( radius.lgCylnOn )
00300 {
00301 fprintf( ioQQQ,
00302 " Luminosity (erg/s) emitted by partial cylinder.\n" );
00303 }
00304
00305 else if( geometry.covgeo == 1. )
00306 {
00307 fprintf( ioQQQ,
00308 " Luminosity (erg/s) emitted by shell with full coverage.\n" );
00309 }
00310
00311 else
00312 {
00313 fprintf( ioQQQ,
00314 " Luminosity (erg/s) emitted by shell with a covering factor of%6.1f%%.\n",
00315 geometry.covgeo*100. );
00316 }
00317 }
00318
00319 else
00320 {
00321 fprintf( ioQQQ, " Intensity (erg/s/cm^2)\n" );
00322 }
00323
00324
00325 fprintf( ioQQQ, " %c\n", ' ' );
00326
00327
00328
00329
00330
00331
00332 if( !atmdat.lgHCaseBOK[1][ipHYDROGEN] )
00333 {
00334 # define NWLH 21
00335
00336 realnum wlh[NWLH] = {6563.e0f, 4861.e0f, 4340.e0f, 4102.e0f, 1.875e4f, 1.282e4f,
00337 1.094e4f, 1.005e4f, 2.625e4f, 2.166e4f, 1.945e4f, 7.458e4f,
00338 4.653e4f, 3.740e4f, 4.051e4f, 7.458e4f, 3.296e4f, 1216.e0f,
00339 1026.e0f, 972.5e0f, 949.7e0f };
00340
00341
00342 for( i=0; i < LineSave.nsum; i++ )
00343 {
00344
00345
00346
00347 if( (strcmp( LineSv[i].chALab,"Ca B" )==0) ||
00348 (strcmp( LineSv[i].chALab,"Ca A" )==0) )
00349 {
00350 realnum errorwave;
00351
00352
00353
00354 errorwave = WavlenErrorGet( LineSv[i].wavelength );
00355 for( j=0; j<NWLH; ++j )
00356 {
00357 if( fabs(LineSv[i].wavelength-wlh[j] ) <= errorwave )
00358 {
00359 LineSv[i].sumlin[0] = 0.;
00360 LineSv[i].sumlin[1] = 0.;
00361 break;
00362 }
00363 }
00364 }
00365 }
00366 }
00367
00368 if( !atmdat.lgHCaseBOK[1][ipHELIUM] )
00369 {
00370
00371 # define NWLHE 20
00372 realnum wlhe[NWLHE] = {1640.e0f, 1215.e0f, 1085.e0f, 1025.e0f, 4686.e0f, 3203.e0f,
00373 2733.e0f, 2511.e0f, 1.012e4f, 6560.e0f, 5412.e0f, 4860.e0f,
00374 1.864e4f, 1.163e4f, 9345.e0f, 8237.e0f, 303.8e0f, 256.3e0f,
00375 243.0e0f, 237.3e0f};
00376 for( i=0; i < LineSave.nsum; i++ )
00377 {
00378 if( (strcmp( LineSv[i].chALab,"Ca B" )==0) ||
00379 (strcmp( LineSv[i].chALab,"Ca A" )==0) )
00380 {
00381 realnum errorwave;
00382
00383
00384
00385 errorwave = WavlenErrorGet( LineSv[i].wavelength );
00386 for( j=0; j<NWLHE; ++j )
00387 {
00388 if( fabs(LineSv[i].wavelength-wlhe[j] ) <= errorwave )
00389 {
00390 LineSv[i].sumlin[0] = 0.;
00391 LineSv[i].sumlin[1] = 0.;
00392 break;
00393 }
00394 }
00395 }
00396 }
00397 }
00398
00399
00400
00401
00402
00403
00404
00405 if( cdLine("TOTL",4861.,&hbeta,&absint)<=0 )
00406 {
00407 if( dense.lgElmtOn[ipHYDROGEN] )
00408 {
00409
00410 fprintf( ioQQQ, " PrtFinal could not find TOTL 4861 with cdLine.\n" );
00411 cdEXIT(EXIT_FAILURE);
00412 }
00413 }
00414
00415 if( dense.lgElmtOn[ipHELIUM] )
00416 {
00417
00418
00419 if( cdLine("He 1",5875.61f,&he5876,&absint)<=0 )
00420 {
00421
00422 if( iso.numLevels_local[ipHE_LIKE][ipHELIUM] >= 14 )
00423 {
00424
00425 fprintf( ioQQQ, " PrtFinal could not find He 1 5876 with cdLine.\n" );
00426 cdEXIT(EXIT_FAILURE);
00427 }
00428 }
00429
00430
00431
00432 if( cdLine("He 2",4686.01f,&he4686,&absint)<=0 )
00433 {
00434
00435 if( iso.numLevels_local[ipH_LIKE][ipHELIUM] > 5 )
00436 {
00437
00438
00439 fprintf( ioQQQ, " PrtFinal could not find He 2 4686 with cdLine.\n" );
00440 cdEXIT(EXIT_FAILURE);
00441 }
00442 }
00443 }
00444 else
00445 {
00446 he5876 = 0.;
00447 absint = 0.;
00448 he4686 = 0.;
00449 }
00450
00451 if( hbeta > 0. )
00452 {
00453 heabun = (he4686*0.078 + he5876*0.739)/hbeta;
00454 }
00455 else
00456 {
00457 heabun = 0.;
00458 }
00459
00460 if( dense.lgElmtOn[ipHELIUM] )
00461 {
00462 hescal = heabun/(dense.gas_phase[ipHELIUM]/dense.gas_phase[ipHYDROGEN]);
00463 }
00464 else
00465 {
00466 hescal = 0.;
00467 }
00468
00469
00470
00471 if( cdLine("O 3",5007.,&o5007,&absint)<=0 )
00472 {
00473 if( dense.lgElmtOn[ipOXYGEN] )
00474 {
00475
00476 fprintf( ioQQQ, " PrtFinal could not find O 3 5007 with cdLine.\n" );
00477 cdEXIT(EXIT_FAILURE);
00478 }
00479 }
00480
00481 if( cdLine("TOTL",4363.,&o4363,&absint)<=0 )
00482 {
00483 if( dense.lgElmtOn[ipOXYGEN] )
00484 {
00485
00486 fprintf( ioQQQ, " PrtFinal could not find TOTL 4363 with cdLine.\n" );
00487 cdEXIT(EXIT_FAILURE);
00488 }
00489 }
00490
00491
00492 if( (o4363 != 0.) && (o5007 != 0.) )
00493 {
00494
00495
00496 bot = o5007 - o4363/oxy.o3enro;
00497
00498 if( bot > 0. )
00499 {
00500 ratio = o4363/bot;
00501 }
00502 else
00503 {
00504 ratio = 0.178;
00505 }
00506
00507 if( ratio > 0.177 )
00508 {
00509
00510 peimbt.toiiir = 1e10;
00511 }
00512 else
00513 {
00514
00515
00516
00517
00518
00519 oxy.o3cs12 = 2.2347f;
00520 oxy.o3cs13 = 0.29811f;
00521 ratio = ratio/1.3333/(oxy.o3cs13/oxy.o3cs12);
00522
00523 ratio = ratio/oxy.o3enro/oxy.o3br32;
00524 peimbt.toiiir = (realnum)(-oxy.o3ex23/log(ratio));
00525 }
00526 }
00527
00528 else
00529 {
00530 peimbt.toiiir = 0.;
00531 }
00532
00533
00534 if( cdLine("Bac ",3646.,&bacobs,&absint)<=0 )
00535 {
00536 fprintf( ioQQQ, " PrtFinal could not find Bac 3546 with cdLine.\n" );
00537 cdEXIT(EXIT_FAILURE);
00538 }
00539
00540
00541 if( hbeta > 0. )
00542 {
00543 bac = bacobs/hbeta;
00544 }
00545 else
00546 {
00547 bac = 0.;
00548 }
00549
00550 if( bac > 0. )
00551 {
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 x = 1.e0/log10(bac);
00571 tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 +
00572 x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));
00573
00574 if( tel > 1. && tel < 5. )
00575 {
00576 peimbt.tbac = (realnum)pow(10.,tel);
00577 }
00578 else
00579 {
00580 peimbt.tbac = 0.;
00581 }
00582 }
00583
00584 else
00585 {
00586 peimbt.tbac = 0.;
00587 }
00588
00589
00590 if( cdLine("thin",3646.,&bacthn,&absint)<=0 )
00591 {
00592 fprintf( ioQQQ, " PrtFinal could not find thin 3646 with cdLine.\n" );
00593 cdEXIT(EXIT_FAILURE);
00594 }
00595
00596
00597 if( cdLine("Ca B",4861.36f,&hbcab,&absint)<=0 )
00598 {
00599 fprintf( ioQQQ, " PrtFinal could not find Ca B 4861 with cdLine.\n" );
00600 cdEXIT(EXIT_FAILURE);
00601 }
00602
00603 if( hbcab > 0. )
00604 {
00605 bacthn /= hbcab;
00606 }
00607 else
00608 {
00609 bacthn = 0.;
00610 }
00611
00612 if( bacthn > 0. )
00613 {
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632 x = 1.e0/log10(bacthn);
00633 tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 +
00634 x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));
00635
00636 if( tel > 1. && tel < 5. )
00637 {
00638 peimbt.tbcthn = (realnum)pow(10.,tel);
00639 }
00640 else
00641 {
00642 peimbt.tbcthn = 0.;
00643 }
00644 }
00645 else
00646 {
00647 peimbt.tbcthn = 0.;
00648 }
00649
00650
00651
00652
00653 peimbt.tohyox = (realnum)((8.5*peimbt.toiiir - 7.5*peimbt.tbcthn - 228200. +
00654 sqrt(POW2(8.5*peimbt.toiiir-7.5*peimbt.tbcthn-228200.)+9.128e5*
00655 peimbt.tbcthn))/2.);
00656
00657 if( peimbt.tohyox > 0. )
00658 {
00659 peimbt.t2hyox = (realnum)((peimbt.tohyox - peimbt.tbcthn)/(1.7*peimbt.tohyox));
00660 }
00661 else
00662 {
00663 peimbt.t2hyox = 0.;
00664 }
00665
00666
00667
00668
00669
00670
00671 scaled = (realnum *)MALLOC( sizeof(realnum)*(unsigned long)LineSave.nsum);
00672
00673
00674 xLog_line_lumin = (double *)MALLOC( sizeof(double)*(unsigned long)LineSave.nsum);
00675
00676
00677
00678 lgPrt = (short int *)MALLOC( sizeof(short int)*(unsigned long)LineSave.nsum);
00679
00680
00681 wavelength = (realnum *)MALLOC( sizeof(realnum)*(unsigned long)LineSave.nsum);
00682
00683
00684 sclsav = (realnum *)MALLOC( sizeof(realnum)*(unsigned long)LineSave.nsum );
00685
00686
00687 chSTyp = (char *)MALLOC( sizeof(char)*(unsigned long)LineSave.nsum );
00688
00689
00690
00691
00692 chSLab = ((char**)MALLOC((unsigned long)LineSave.nsum*sizeof(char*)));
00693
00694
00695 for( i=0; i<LineSave.nsum; ++i)
00696 {
00697 chSLab[i] = (char*)MALLOC(5*sizeof(char));
00698 }
00699
00700
00701 ipSortLines = (long *)MALLOC( sizeof(long)*(unsigned long)LineSave.nsum );
00702
00703 ASSERT( LineSave.ipNormWavL >= 0 );
00704 for( ipEmType=0; ipEmType<2; ++ipEmType )
00705 {
00706
00707
00708 snorm = LineSv[LineSave.ipNormWavL].sumlin[ipEmType];
00709
00710
00711 if( ((snorm <= SMALLDOUBLE ) || (LineSave.ipNormWavL < 0)) || (LineSave.ipNormWavL > LineSave.nsum) )
00712 {
00713 fprintf( ioQQQ, "\n\n >>PROBLEM Normalization line has small or zero intensity, its value was %.2e and its intensity was set to 1."
00714 "\n >>Please consider using another normalization line (this is set with the NORMALIZE command).\n" , snorm);
00715 fprintf( ioQQQ, " >>The relative intensities will be meaningless, and many lines may appear too faint.\n" );
00716 snorm = 1.;
00717 }
00718 for( i=0; i < LineSave.nsum; i++ )
00719 {
00720
00721
00722
00723
00724
00725 double scale = LineSv[i].sumlin[ipEmType]/snorm*LineSave.ScaleNormLine;
00726
00727 scale = MIN2(BIGFLOAT , scale );
00728 scale = MAX2( -BIGFLOAT , scale );
00729
00730
00731 scaled[i] = (realnum)scale;
00732
00733 if( LineSv[i].sumlin[ipEmType] > 0. )
00734 {
00735 xLog_line_lumin[i] = log10(LineSv[i].sumlin[ipEmType]) + radius.Conv2PrtInten;
00736 }
00737 else
00738 {
00739 xLog_line_lumin[i] = -38.;
00740 }
00741 }
00742
00743
00744 for( i=0; i < LineSave.nsum; i++ )
00745 {
00746
00747 if( strcmp(LineSv[i].chALab,"Unit") == 0 )
00748 {
00749 lgPrt[i] = false;
00750 }
00751 else if( strcmp(LineSv[i].chALab,"Coll") == 0 && !prt.lgPrnColl )
00752 {
00753 lgPrt[i] = false;
00754 }
00755 else if( strcmp(LineSv[i].chALab,"Pump") == 0 && !prt.lgPrnPump )
00756 {
00757 lgPrt[i] = false;
00758 }
00759 else if( strncmp(LineSv[i].chALab,"Inw",3) == 0 && !prt.lgPrnInwd )
00760 {
00761 lgPrt[i] = false;
00762 }
00763 else if( strcmp(LineSv[i].chALab,"Heat") == 0 && !prt.lgPrnHeat )
00764 {
00765 lgPrt[i] = false;
00766 }
00767
00768
00769 else if( !prt.lgPrnDiff &&
00770 ( (strcmp(LineSv[i].chALab,"nFnu") == 0) || (strcmp(LineSv[i].chALab,"nInu") == 0) ) )
00771 {
00772 lgPrt[i] = false;
00773 }
00774 else
00775 {
00776 lgPrt[i] = true;
00777 }
00778 }
00779
00780
00781 nNegIntenLines = 0;
00782
00783
00784 for(i=0; i< 32; i++ )
00785 {
00786 ipNegIntensity[i] = LONG_MAX;
00787 }
00788
00789 for(i=0;i<8;++i)
00790 {
00791 d[i] = -DBL_MAX;
00792 }
00793
00794
00795 const char chIntensityType[2][10]={"Intrinsic" , " Emergent"};
00796 ASSERT( ipEmType==0 || ipEmType==1 );
00797
00798
00799 fprintf( ioQQQ, "\n" );
00800 if( prt.lgPrtLineArray )
00801 fprintf( ioQQQ, " " );
00802 fprintf( ioQQQ, "%s" , chIntensityType[ipEmType] );
00803 fprintf( ioQQQ, " line intensities\n" );
00804
00805
00806 if( ipEmType==1 && iteration==1 )
00807 fprintf(ioQQQ," Caution: emergent intensities are not reliable on the "
00808 "first iteration.\n");
00809
00810
00811 if( prt.lgFaintOn )
00812 {
00813 iprnt = 0;
00814 for( i=0; i < LineSave.nsum; i++ )
00815 {
00816
00817 ASSERT( iprnt <= i);
00818 if( scaled[i] >= prt.TooFaint && lgPrt[i] )
00819 {
00820 if( prt.lgPrtLineLog )
00821 {
00822 xLog_line_lumin[iprnt] = log10(LineSv[i].sumlin[ipEmType]) + radius.Conv2PrtInten;
00823 }
00824 else
00825 {
00826 xLog_line_lumin[iprnt] = LineSv[i].sumlin[ipEmType] * pow(10.,radius.Conv2PrtInten);
00827 }
00828 sclsav[iprnt] = scaled[i];
00829 chSTyp[iprnt] = LineSv[i].chSumTyp;
00830
00831
00832 ASSERT( strlen( LineSv[i].chALab )<5 );
00833 strcpy( chSLab[iprnt], LineSv[i].chALab );
00834 wavelength[iprnt] = LineSv[i].wavelength;
00835 ++iprnt;
00836 }
00837 else if( -scaled[i] > prt.TooFaint && lgPrt[i] )
00838 {
00839
00840 ipNegIntensity[nNegIntenLines] = i;
00841 nNegIntenLines = MIN2(32,nNegIntenLines+1);
00842 }
00843
00844
00845 else if( strcmp( LineSv[i].chALab, "####" ) == 0 &&!prt.lgSortLines )
00846 {
00847 strcpy( chSLab[iprnt], LineSv[i].chALab );
00848 xLog_line_lumin[iprnt] = 0.;
00849 sclsav[iprnt] = 0.;
00850 chSTyp[iprnt] = LineSv[i].chSumTyp;
00851 wavelength[iprnt] = LineSv[i].wavelength;
00852 ++iprnt;
00853 }
00854 }
00855 }
00856
00857 else
00858 {
00859
00860 iprnt = LineSave.nsum;
00861 for( i=0; i < LineSave.nsum; i++ )
00862 {
00863 if( strcmp( LineSv[i].chALab, "####" ) == 0 )
00864 {
00865 strcpy( chSLab[i], LineSv[i].chALab );
00866 xLog_line_lumin[i] = 0.;
00867 sclsav[i] = 0.;
00868 chSTyp[i] = LineSv[i].chSumTyp;
00869 wavelength[i] = LineSv[i].wavelength;
00870 }
00871 else
00872 {
00873 sclsav[i] = scaled[i];
00874 chSTyp[i] = LineSv[i].chSumTyp;
00875 strcpy( chSLab[i], LineSv[i].chALab );
00876 wavelength[i] = LineSv[i].wavelength;
00877 }
00878 if( scaled[i] < 0. )
00879 {
00880 ipNegIntensity[nNegIntenLines] = i;
00881 nNegIntenLines = MIN2(32,nNegIntenLines+1);
00882 }
00883 }
00884 }
00885
00886
00887 ASSERT( iprnt > 0. );
00888
00889
00890
00891 if( prt.lgSortLines )
00892 {
00893 int lgFlag;
00894 if( prt.lgSortLineWavelength )
00895 {
00896
00897 if( prt.wlSort1 >-0.1 )
00898 {
00899 j = 0;
00900
00901 for( i=0; i<iprnt; ++i )
00902 {
00903 if( wavelength[i]>= prt.wlSort1 && wavelength[i]<= prt.wlSort2 )
00904 {
00905 if( j!=i )
00906 {
00907 sclsav[j] = sclsav[i];
00908 chSTyp[j] = chSTyp[i];
00909 strcpy( chSLab[j], chSLab[i] );
00910 wavelength[j] = wavelength[i];
00911
00912
00913 xLog_line_lumin[j] = xLog_line_lumin[i];
00914 }
00915 ++j;
00916 }
00917 }
00918 iprnt = j;
00919 }
00920
00921 spsort(wavelength,
00922 iprnt,
00923 ipSortLines,
00924
00925
00926 1,
00927 &lgFlag);
00928 if( lgFlag )
00929 TotalInsanity();
00930 }
00931 else if( prt.lgSortLineIntensity )
00932 {
00933
00934 spsort(sclsav,
00935 iprnt,
00936 ipSortLines,
00937
00938
00939 -1,
00940 &lgFlag);
00941 if( lgFlag )
00942 TotalInsanity();
00943 }
00944 else
00945 {
00946
00947 TotalInsanity();
00948 }
00949 }
00950 else
00951 {
00952
00953 for( i=0; i<iprnt; ++i )
00954 {
00955 ipSortLines[i] = i;
00956 }
00957 }
00958
00959
00960
00961
00962
00963 if( prt.lgPrtLineArray )
00964 {
00965
00966 if( LineSave.sig_figs >= 5 )
00967 {
00968 nline = (iprnt + 2)/3;
00969 }
00970 else
00971 {
00972
00973
00974 nline = (iprnt + 3)/4;
00975 }
00976 }
00977 else
00978 {
00979
00980 nline = iprnt;
00981 }
00982
00983
00984 for( i=0; i < nline; i++ )
00985 {
00986 fprintf( ioQQQ, " " );
00987
00988
00989
00990 for( j=i; j<iprnt; j = j + nline)
00991 {
00992
00993 long ipLin = ipSortLines[j];
00994
00995
00996 if( strcmp( chSLab[ipLin], "####" ) == 0 )
00997 {
00998
00999
01000
01001 fprintf( ioQQQ, "%s",LineSave.chHoldComments[(int)wavelength[ipLin]] );
01002 }
01003 else
01004 {
01005
01006 fprintf( ioQQQ, "%4.4s ",chSLab[ipLin] );
01007
01008
01009 prt_wl( ioQQQ , wavelength[ipLin]);
01010
01011
01012
01013 if( prt.lgPrtLineLog )
01014 {
01015 fprintf( ioQQQ, " %7.3f", xLog_line_lumin[ipLin] );
01016 }
01017 else
01018 {
01019
01020 fprintf( ioQQQ, " %.2e ", xLog_line_lumin[ipLin] );
01021 }
01022
01023
01024
01025
01026 if( sclsav[ipLin] < 9999.5 )
01027 {
01028 fprintf( ioQQQ, "%9.4f",sclsav[ipLin] );
01029 }
01030 else if( sclsav[ipLin] < 99999.5 )
01031 {
01032 fprintf( ioQQQ, "%9.3f",sclsav[ipLin] );
01033 }
01034 else if( sclsav[ipLin] < 999999.5 )
01035 {
01036 fprintf( ioQQQ, "%9.2f",sclsav[ipLin] );
01037 }
01038 else if( sclsav[ipLin] < 9999999.5 )
01039 {
01040 fprintf( ioQQQ, "%9.1f",sclsav[ipLin] );
01041 }
01042 else
01043 {
01044 fprintf( ioQQQ, "*********" );
01045 }
01046
01047
01048 fprintf( ioQQQ, " " );
01049 }
01050 }
01051
01052 fprintf( ioQQQ, " \n" );
01053 }
01054
01055 if( nNegIntenLines > 0 )
01056 {
01057 fprintf( ioQQQ, " Lines with negative intensities - Linear intensities relative to normalize line.\n" );
01058 fprintf( ioQQQ, " " );
01059
01060 for( i=0; i < nNegIntenLines; ++i )
01061 {
01062 fprintf( ioQQQ, "%ld %s %.0f %.1e, ",
01063 ipNegIntensity[i],
01064 LineSv[ipNegIntensity[i]].chALab,
01065 LineSv[ipNegIntensity[i]].wavelength,
01066 scaled[ipNegIntensity[i]] );
01067 }
01068 fprintf( ioQQQ, "\n" );
01069 }
01070 }
01071
01072
01073 j = 0;
01074 for( i=0; i < LineSave.nsum; i++ )
01075 {
01076 a = (double)LineSv[i].sumlin[LineSave.lgLineEmergent]/(double)thermal.totcol;
01077 if( (a >= 0.05) && LineSv[i].chSumTyp == 'c' )
01078 {
01079 ipNegIntensity[j] = i;
01080 d[j] = a;
01081 j = MIN2(j+1,7);
01082 }
01083 }
01084
01085 fprintf( ioQQQ, "\n\n\n %s\n", input.chTitle );
01086 if( j != 0 )
01087 {
01088 fprintf( ioQQQ, " Cooling: " );
01089 for( i=0; i < j; i++ )
01090 {
01091 fprintf( ioQQQ, " %4.4s ",
01092 LineSv[ipNegIntensity[i]].chALab);
01093
01094 prt_wl( ioQQQ, LineSv[ipNegIntensity[i]].wavelength );
01095
01096 fprintf( ioQQQ, ":%5.3f",
01097 d[i] );
01098 }
01099 fprintf( ioQQQ, " \n" );
01100 }
01101
01102
01103 j = 0;
01104 for( i=0; i < LineSave.nsum; i++ )
01105 {
01106 a = (double)LineSv[i].sumlin[LineSave.lgLineEmergent]/(double)thermal.power;
01107 if( (a >= 0.05) && LineSv[i].chSumTyp == 'h' )
01108 {
01109 ipNegIntensity[j] = i;
01110 d[j] = a;
01111 j = MIN2(j+1,7);
01112 }
01113 }
01114
01115 if( j != 0 )
01116 {
01117 fprintf( ioQQQ, " Heating: " );
01118 for( i=0; i < j; i++ )
01119 {
01120 fprintf( ioQQQ, " %4.4s ",
01121 LineSv[ipNegIntensity[i]].chALab);
01122
01123 prt_wl(ioQQQ, LineSv[ipNegIntensity[i]].wavelength);
01124
01125 fprintf( ioQQQ, ":%5.3f",
01126 d[i] );
01127 }
01128 fprintf( ioQQQ, " \n" );
01129 }
01130
01131
01132 qlow = 0.;
01133 plow = 0.;
01134 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
01135 {
01136
01137 plow += (rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i])*
01138 EN1RYD*rfield.anu[i];
01139 qlow += rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i];
01140 }
01141
01142 qlow *= radius.r1r0sq;
01143 plow *= radius.r1r0sq;
01144 if( plow > 0. )
01145 {
01146 qlow = log10(qlow) + radius.Conv2PrtInten;
01147 plow = log10(plow) + radius.Conv2PrtInten;
01148 }
01149 else
01150 {
01151 qlow = 0.;
01152 plow = 0.;
01153 }
01154
01155 qion = 0.;
01156 pion = 0.;
01157 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ )
01158 {
01159
01160 pion += (rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i])*
01161 EN1RYD*rfield.anu[i];
01162 qion += rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i];
01163 }
01164
01165 qion *= radius.r1r0sq;
01166 pion *= radius.r1r0sq;
01167
01168 if( pion > 0. )
01169 {
01170 qion = log10(qion) + radius.Conv2PrtInten;
01171 pion = log10(pion) + radius.Conv2PrtInten;
01172 }
01173 else
01174 {
01175 qion = 0.;
01176 pion = 0.;
01177 }
01178
01179
01180 if( rfield.qhtot > 0. )
01181 {
01182 if( rfield.lgUSphON )
01183 {
01184
01185 usp = rfield.qhtot/POW2(rfield.rstrom/radius.rinner)/
01186 2.998e10/dense.gas_phase[ipHYDROGEN];
01187 usp = log10(usp);
01188 }
01189 else
01190 {
01191
01192 usp = rfield.qhtot/radius.r1r0sq/2.998e10/dense.gas_phase[ipHYDROGEN];
01193 usp = log10(usp);
01194 }
01195 }
01196
01197 else
01198 {
01199 usp = -37.;
01200 }
01201
01202 if( rfield.uh > 0. )
01203 {
01204 uhl = log10(rfield.uh);
01205 }
01206 else
01207 {
01208 uhl = -37.;
01209 }
01210
01211 if( rfield.uheii > 0. )
01212 {
01213 uhel = log10(rfield.uheii);
01214 }
01215 else
01216 {
01217 uhel = -37.;
01218 }
01219
01220 fprintf( ioQQQ,
01221 "\n IONIZE PARMET: U(1-)%8.4f U(4-):%8.4f U(sp):%6.2f "
01222 "Q(ion): %7.3f L(ion)%7.3f Q(low):%7.3f L(low)%7.3f\n",
01223 uhl, uhel, usp, qion, pion, qlow, plow );
01224
01225 a = fabs((thermal.power-thermal.totcol)*100.)/thermal.power;
01226
01227 if( thermal.power > 0. )
01228 {
01229 powerl = log10(thermal.power) + radius.Conv2PrtInten;
01230 }
01231 else
01232 {
01233 powerl = 0.;
01234 }
01235
01236 if( thermal.totcol > 0. )
01237 {
01238 thermal.totcol = log10(thermal.totcol) + radius.Conv2PrtInten;
01239 }
01240 else
01241 {
01242 thermal.totcol = 0.;
01243 }
01244
01245 if( thermal.FreeFreeTotHeat > 0. )
01246 {
01247 thermal.FreeFreeTotHeat = log10(thermal.FreeFreeTotHeat) + radius.Conv2PrtInten;
01248 }
01249 else
01250 {
01251 thermal.FreeFreeTotHeat = 0.;
01252 }
01253
01254
01255 reclin = totlin('r');
01256 if( reclin > 0. )
01257 {
01258 reclin = log10(reclin) + radius.Conv2PrtInten;
01259 }
01260 else
01261 {
01262 reclin = 0.;
01263 }
01264
01265 fprintf( ioQQQ,
01266 " ENERGY BUDGET: Heat:%8.3f Coolg:%8.3f Error:%5.1f%% Rec Lin:%8.3f F-F H%7.3f P(rad/tot)max ",
01267 powerl,
01268 thermal.totcol,
01269 a,
01270 reclin,
01271 thermal.FreeFreeTotHeat );
01272 PrintE82( ioQQQ, pressure.RadBetaMax );
01273 fprintf( ioQQQ, "\n");
01274
01275
01276
01277 coleff = opac.TauAbsGeo[0][rt.ipxry-1] - rt.tauxry;
01278 coleff /= 2.14e-22;
01279
01280
01281 gett2(&peimbt.t2hstr);
01282
01283
01284 gett2o3(&peimbt.t2o3str);
01285
01286 fprintf( ioQQQ, "\n Col(Heff): ");
01287 PrintE93(ioQQQ, coleff);
01288 fprintf( ioQQQ, " snd travl time ");
01289 PrintE82(ioQQQ, timesc.sound);
01290 fprintf( ioQQQ, " sec Te-low: ");
01291 PrintE82(ioQQQ, thermal.tlowst);
01292 fprintf( ioQQQ, " Te-hi: ");
01293 PrintE82(ioQQQ, thermal.thist);
01294
01295
01296
01297 fprintf( ioQQQ, " G0TH85: ");
01298 PrintE82( ioQQQ, hmi.UV_Cont_rel2_Habing_TH85_face );
01299
01300
01301 fprintf( ioQQQ, " G0DB96:");
01302 PrintE82( ioQQQ, hmi.UV_Cont_rel2_Draine_DB96_face );
01303
01304 fprintf( ioQQQ, "\n");
01305
01306 fprintf( ioQQQ, " Emiss Measure n(e)n(p) dl ");
01307 PrintE93(ioQQQ, colden.dlnenp);
01308 fprintf( ioQQQ, " n(e)n(He+)dl ");
01309 PrintE93(ioQQQ, colden.dlnenHep);
01310 fprintf( ioQQQ, " En(e)n(He++) dl ");
01311 PrintE93(ioQQQ, colden.dlnenHepp);
01312 fprintf( ioQQQ, " ne nC+:");
01313 PrintE82(ioQQQ, colden.dlnenCp);
01314 fprintf( ioQQQ, "\n");
01315
01316
01317 if( rfield.EnergyBremsThin > 0. )
01318 {
01319 bremtk = RYDLAM*1e-8/rfield.EnergyBremsThin;
01320 }
01321 else
01322 {
01323 bremtk = 1e30;
01324 }
01325
01326
01327 fprintf( ioQQQ, " He/Ha:");
01328 PrintE82( ioQQQ, heabun);
01329
01330
01331 fprintf(ioQQQ, " =%7.2f*true Lthin:",hescal);
01332
01333
01334 PrintE82( ioQQQ, bremtk);
01335
01336
01337
01338
01339
01340
01341 if( nzone > 0 )
01342 {
01343
01344
01345
01346
01347 znit = (double)(conv.nTotalIoniz-conv.nTotalIoniz_start)/(double)(nzone);
01348 }
01349 else
01350 {
01351 znit = 0.;
01352 }
01353
01354 fprintf(ioQQQ, " itr/zn:%5.2f",znit);
01355
01356
01357 fprintf(ioQQQ, " H2 itr/zn:%6.2f",H2_itrzn());
01358
01359
01360 fprintf(ioQQQ, " File Opacity: F" );
01361
01362
01363 {
01364
01365 double xmass = log10( SDIV(dense.xMassTotal) );
01366 xmass += (realnum)(1.0992099 + 2.*log10(radius.rinner));
01367 fprintf(ioQQQ," Mass tot %.3f",
01368 xmass);
01369 }
01370 fprintf(ioQQQ,"\n");
01371
01372
01373
01374
01375 if( cdTemp( "opti",0,&THI,"volume" ) )
01376 {
01377 fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm opt.\n");
01378 TotalInsanity();
01379 }
01380 fprintf( ioQQQ, " Temps(21 cm) T(21cm/Ly a) ");
01381 PrintE82(ioQQQ, THI );
01382
01383
01384
01385
01386
01387 if( cdTemp( "21cm",0,&THI,"radius" ) )
01388 {
01389 fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm temp.\n");
01390 TotalInsanity();
01391 }
01392 fprintf(ioQQQ, " T(<nH/Tkin>) ");
01393 PrintE82(ioQQQ, THI);
01394
01395
01396
01397 if( cdTemp( "spin",0,&THI,"radius" ) )
01398 {
01399 fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm spin.\n");
01400 TotalInsanity();
01401 }
01402 fprintf(ioQQQ, " T(<nH/Tspin>) ");
01403 PrintE82(ioQQQ, THI);
01404
01405
01406 THI *= (1. - sexp( HFLines[0].Emis->TauIn ) );
01407 fprintf( ioQQQ, " TB21cm:");
01408 PrintE82(ioQQQ, THI);
01409
01410
01411 if( cdTemp( "spin",0,&THI,"volume" ) )
01412 {
01413 fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting spin temp.\n");
01414 TotalInsanity();
01415 }
01416 fprintf( ioQQQ, "\n <Tspin> ");
01417 PrintE82(ioQQQ, THI);
01418 fprintf( ioQQQ, " N(H0/Tspin) ");
01419 PrintE82(ioQQQ, colden.H0_ov_Tspin );
01420 fprintf( ioQQQ, " N(OH/Tkin) ");
01421 PrintE82(ioQQQ, colden.OH_ov_Tspin );
01422
01423
01424 bot = cdB21cm();
01425 fprintf( ioQQQ, " B(21cm):");
01426 PrintE82(ioQQQ, bot );
01427
01428
01429 fprintf(ioQQQ, "\n");
01430
01431
01432 rate = timesc.TimeErode*2e-26;
01433 if( rate > SMALLFLOAT )
01434 {
01435 ferode = 1./rate;
01436 }
01437 else
01438 {
01439 ferode = 0.;
01440 }
01441
01442
01443 if( wind.acldr > 0. )
01444 {
01445 wind.AccelAver /= wind.acldr;
01446 }
01447 else
01448 {
01449 wind.AccelAver = 0.;
01450 }
01451
01452
01453 wmean = colden.wmas/SDIV(colden.TotMassColl);
01454
01455 dmean = colden.TotMassColl/SDIV(radius.depth_x_fillfac);
01456 tmean = colden.tmas/SDIV(colden.TotMassColl);
01457
01458 wmean = colden.wmas/SDIV(colden.TotMassColl);
01459
01460 fprintf( ioQQQ, " <a>:");
01461 PrintE82(ioQQQ , wind.AccelAver);
01462 fprintf( ioQQQ, " erdeFe");
01463 PrintE71(ioQQQ , ferode);
01464 fprintf( ioQQQ, " Tcompt");
01465 PrintE82(ioQQQ , timesc.tcmptn);
01466 fprintf( ioQQQ, " Tthr");
01467 PrintE82(ioQQQ , timesc.time_therm_long);
01468 fprintf( ioQQQ, " <Tden>: ");
01469 PrintE82(ioQQQ , tmean);
01470 fprintf( ioQQQ, " <dens>:");
01471 PrintE82(ioQQQ , dmean);
01472 fprintf( ioQQQ, " <MolWgt>");
01473 PrintE82(ioQQQ , wmean);
01474 fprintf(ioQQQ,"\n");
01475
01476
01477 if( tmean > 0. )
01478 {
01479 rjeans = 7.79637 + (log10(tmean) - log10(dense.wmole) - log10(dense.xMassDensity*
01480 geometry.FillFac))/2.;
01481 }
01482 else
01483 {
01484
01485 rjeans = 40.f;
01486 }
01487
01488 if( rjeans < 36. )
01489 {
01490 rjeans = (double)pow(10.,rjeans);
01491
01492 ajmass = 3.*log10(rjeans/2.) + log10(4.*PI/3.*dense.xMassDensity*
01493 geometry.FillFac) - log10(SOLAR_MASS);
01494 if( ajmass < 37 )
01495 {
01496 ajmass = pow(10.,ajmass);
01497 }
01498 else
01499 {
01500 ajmass = 0.;
01501 }
01502 }
01503 else
01504 {
01505 rjeans = 0.;
01506 ajmass = 0.;
01507 }
01508
01509
01510 ajmin = colden.ajmmin - log10(SOLAR_MASS);
01511 if( ajmin < 37 )
01512 {
01513 ajmin = pow(10.,ajmin);
01514 }
01515 else
01516 {
01517 ajmin = 0.;
01518 }
01519
01520
01521 if( rfield.anu[rfield.nflux-1] > 150. )
01522 {
01523
01524 ip2kev = ipoint(147.);
01525 ip2500 = ipoint(0.3648);
01526
01527
01528 bottom = rfield.flux[ip2500-1]*
01529 rfield.anu[ip2500-1]/rfield.widflx[ip2500-1];
01530
01531 top = rfield.flux[ip2kev-1]*
01532 rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1];
01533
01534
01535 if( bottom > 1e-30 && top > 1e-30 )
01536 {
01537 ratio = log10(top) - log10(bottom);
01538 if( ratio < 36. && ratio > -36. )
01539 {
01540 ratio = pow(10.,ratio);
01541 }
01542 else
01543 {
01544 ratio = 0.;
01545 }
01546 }
01547
01548 else
01549 {
01550 ratio = 0.;
01551 }
01552
01553 if( ratio > 0. )
01554 {
01555
01556 double freq_ratio = rfield.anu[ip2kev-1]/rfield.anu[ip2500-1];
01557 alfox = log(ratio)/log(freq_ratio);
01558 }
01559 else
01560 {
01561 alfox = 0.;
01562 }
01563 }
01564 else
01565 {
01566 alfox = 0.;
01567 }
01568
01569 if( colden.rjnmin < 37 )
01570 {
01571 colden.rjnmin = (realnum)pow((realnum)10.f,colden.rjnmin);
01572 }
01573 else
01574 {
01575 colden.rjnmin = FLT_MAX;
01576 }
01577
01578
01579
01580 fprintf( ioQQQ, " Mean Jeans l(cm)");
01581 PrintE82(ioQQQ, rjeans);
01582 fprintf( ioQQQ, " M(sun)");
01583 PrintE82(ioQQQ, ajmass);
01584 fprintf( ioQQQ, " smallest: len(cm):");
01585 PrintE82(ioQQQ, colden.rjnmin);
01586 fprintf( ioQQQ, " M(sun):");
01587 PrintE82(ioQQQ, ajmin);
01588 fprintf( ioQQQ, " a_ox tran:%6.2f\n" , alfox);
01589
01590 if( prt.lgPrintTime )
01591 {
01592
01593 alfox = cdExecTime();
01594 }
01595 else
01596 {
01597
01598
01599 alfox = 0.;
01600 }
01601
01602
01603 fprintf( ioQQQ,
01604 " Hatom level%3ld HatomType:%4.4s HInducImp %2c"
01605 " He tot level:%3ld He2 level: %3ld ExecTime%8.1f\n",
01606
01607 iso.numLevels_local[ipH_LIKE][ipHYDROGEN],
01608 hydro.chHTopType,
01609 TorF(hydro.lgHInducImp),
01610
01611 iso.numLevels_local[ipHE_LIKE][ipHELIUM],
01612
01613 iso.numLevels_local[ipH_LIKE][ipHELIUM] ,
01614 alfox );
01615
01616
01617 fprintf( ioQQQ,
01618 " ConvrgError(%%) <eden>%7.3f MaxEden%7.3f <H-C>%7.2f Max(H-C)%8.2f <Press>%8.3f MaxPrs er%7.3f\n",
01619 conv.AverEdenError/nzone*100. ,
01620 conv.BigEdenError*100. ,
01621 conv.AverHeatCoolError/nzone*100. ,
01622 conv.BigHeatCoolError*100. ,
01623 conv.AverPressError/nzone*100. ,
01624 conv.BigPressError*100. );
01625
01626 fprintf(ioQQQ,
01627 " Continuity(%%) chng Te%6.1f elec den%6.1f n(H2)%7.1f n(CO) %7.1f\n",
01628 phycon.BigJumpTe*100. ,
01629 phycon.BigJumpne*100. ,
01630 phycon.BigJumpH2*100. ,
01631 phycon.BigJumpCO*100. );
01632
01633
01634 aver("prin",0.,0.," ");
01635
01636
01637
01638 if( dense.gas_phase[ipHYDROGEN] < peimbt.tsqden )
01639 {
01640 fprintf( ioQQQ, " \n" );
01641
01642
01643 fprintf( ioQQQ, " Peimbert T(OIIIr)");
01644 PrintE82( ioQQQ , peimbt.toiiir);
01645
01646
01647 fprintf( ioQQQ, " T(Bac)");
01648 PrintE82( ioQQQ , peimbt.tbac);
01649
01650
01651 fprintf( ioQQQ, " T(Hth)");
01652 PrintE82( ioQQQ , peimbt.tbcthn);
01653
01654
01655 fprintf( ioQQQ, " t2(Hstrc)");
01656 fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hstr));
01657
01658
01659 fprintf( ioQQQ, " T(O3-BAC)");
01660 PrintE82( ioQQQ , peimbt.tohyox);
01661
01662
01663 fprintf( ioQQQ, " t2(O3-BC)");
01664 fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hyox));
01665
01666
01667 fprintf( ioQQQ, " t2(O3str)");
01668 fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2o3str));
01669
01670 fprintf( ioQQQ, "\n");
01671
01672 if( gv.lgDustOn )
01673 {
01674 fprintf( ioQQQ, " Be careful: grains exist. This spectrum was not corrected for reddening before analysis.\n" );
01675 }
01676 }
01677
01678
01679
01680 if( gv.lgDustOn && gv.lgGrainPhysicsOn )
01681 {
01682 long int i0,
01683 i1;
01684 char chQHeat;
01685 double AV , AB;
01686 double total_dust2gas = 0.;
01687
01688 fprintf( ioQQQ, "\n Average Grain Properties (over radius):\n" );
01689
01690 for( i0=0; i0 < gv.nBin; i0 += 10 )
01691 {
01692 if( i0 > 0 )
01693 fprintf( ioQQQ, "\n" );
01694
01695
01696 i1 = MIN2(i0+10,gv.nBin);
01697
01698 fprintf( ioQQQ, " " );
01699 for( nd=i0; nd < i1; nd++ )
01700 {
01701 chQHeat = (char)(gv.bin[nd]->lgEverQHeat ? '*' : ' ');
01702 fprintf( ioQQQ, "%-12.12s%c", gv.bin[nd]->chDstLab, chQHeat );
01703 }
01704 fprintf( ioQQQ, "\n" );
01705
01706 fprintf( ioQQQ, " nd:" );
01707 for( nd=i0; nd < i1; nd++ )
01708 {
01709 if( nd != i0 ) fprintf( ioQQQ," " );
01710 fprintf( ioQQQ, "%7ld ", nd );
01711 }
01712 fprintf( ioQQQ, "\n" );
01713
01714 fprintf( ioQQQ, " <Tgr>:" );
01715 for( nd=i0; nd < i1; nd++ )
01716 {
01717 if( nd != i0 ) fprintf( ioQQQ," " );
01718 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdust/radius.depth_x_fillfac));
01719 }
01720 fprintf( ioQQQ, "\n" );
01721
01722 fprintf( ioQQQ, " <Vel>:" );
01723 for( nd=i0; nd < i1; nd++ )
01724 {
01725 if( nd != i0 ) fprintf( ioQQQ," " );
01726 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdft/radius.depth_x_fillfac));
01727 }
01728 fprintf( ioQQQ, "\n" );
01729
01730 fprintf( ioQQQ, " <Pot>:" );
01731 for( nd=i0; nd < i1; nd++ )
01732 {
01733 if( nd != i0 ) fprintf( ioQQQ," " );
01734 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdpot/radius.depth_x_fillfac));
01735 }
01736 fprintf( ioQQQ, "\n" );
01737
01738 fprintf( ioQQQ, " <D/G>:" );
01739 for( nd=i0; nd < i1; nd++ )
01740 {
01741 if( nd != i0 ) fprintf( ioQQQ," " );
01742 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avDGRatio/radius.depth_x_fillfac));
01743
01744 total_dust2gas += gv.bin[nd]->avDGRatio/radius.depth_x_fillfac;
01745 }
01746 fprintf( ioQQQ, "\n" );
01747 }
01748
01749 fprintf(ioQQQ," Dust to gas ratio (by mass):");
01750 fprintf(ioQQQ,PrintEfmt("%10.3e", total_dust2gas));
01751
01752
01753
01754
01755
01756 AV = rfield.extin_mag_V_point/SDIV(colden.colden[ipCOL_HTOT]);
01757 AB = rfield.extin_mag_B_point/SDIV(colden.colden[ipCOL_HTOT]);
01758
01759 fprintf(ioQQQ,", A(V)/N(H)(pnt):%.3e, (ext):%.3e",
01760 AV,
01761 rfield.extin_mag_V_extended/SDIV(colden.colden[ipCOL_HTOT]) );
01762
01763
01764 fprintf(ioQQQ,", R:");
01765
01766
01767 if( fabs(AB-AV)>SMALLFLOAT )
01768 {
01769 fprintf(ioQQQ,"%.3e", AV/(AB-AV) );
01770 }
01771 else
01772 {
01773 fprintf(ioQQQ,"%.3e", 0. );
01774 }
01775 fprintf(ioQQQ," AV(ext):%.3e (pnt):%.3e\n",
01776 rfield.extin_mag_V_extended,
01777 rfield.extin_mag_V_point);
01778 }
01779
01780
01781 free( wavelength );
01782 free( ipSortLines );
01783 free( sclsav );
01784 free( lgPrt );
01785 free( chSTyp );
01786
01787
01788 for( i=0; i < LineSave.nsum; ++i )
01789 {
01790 free( chSLab[i] );
01791 }
01792 free( chSLab );
01793
01794 free( scaled );
01795 free( xLog_line_lumin );
01796
01797
01798 if( !prt.lgPrtShort && called.lgTalk )
01799 {
01800
01801
01802 PrtAllTau();
01803
01804
01805 if( iterations.lgLastIt )
01806 {
01807
01808 if( prt.lgPrintColumns )
01809 PrtColumns(ioQQQ);
01810
01811 PrtMeanIon('i', false, ioQQQ);
01812
01813 PrtMeanIon('i', true , ioQQQ);
01814
01815 PrtMeanIon('t', false , ioQQQ);
01816
01817 PrtMeanIon('t', true , ioQQQ);
01818
01819 PrtContinuum();
01820 }
01821 }
01822
01823
01824 fprintf( ioQQQ, "%s\n\n", input.chTitle );
01825 return;
01826 }
01827
01828
01829
01830 long int StuffComment( const char * chComment )
01831 {
01832 long int n , i;
01833
01834 DEBUG_ENTRY( "StuffComment()" );
01835
01836
01837 if( LineSave.ipass == 0 )
01838 {
01839 if( LineSave.nComment >= NHOLDCOMMENTS )
01840 {
01841 fprintf( ioQQQ, " Too many comments have been entered into line array; increase the value of NHOLDCOMMENTS.\n" );
01842 cdEXIT(EXIT_FAILURE);
01843 }
01844
01845
01846 # define NWIDTH 33
01847 strcpy( LineSave.chHoldComments[LineSave.nComment], chComment );
01848
01849
01850 n = (long)strlen( LineSave.chHoldComments[LineSave.nComment] );
01851 for( i=0; i<NWIDTH-n-1-6; ++i )
01852 {
01853 strcat( LineSave.chHoldComments[LineSave.nComment], ".");
01854 }
01855
01856 strcat( LineSave.chHoldComments[LineSave.nComment], "..");
01857
01858 for( i=0; i<6; ++i )
01859 {
01860 strcat( LineSave.chHoldComments[LineSave.nComment], " ");
01861 }
01862 }
01863
01864 ++LineSave.nComment;
01865 return( LineSave.nComment-1 );
01866 }
01867
01868
01869 STATIC void gett2(realnum *tsqr)
01870 {
01871 long int i;
01872
01873 double tmean;
01874 double a,
01875 as,
01876 b;
01877
01878 DEBUG_ENTRY( "gett2()" );
01879
01880
01881 a = 0.;
01882 b = 0.;
01883
01884 ASSERT( nzone < struc.nzlim);
01885 for( i=0; i < nzone; i++ )
01886 {
01887 as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
01888 (double)(struc.ednstr[i]);
01889 a += (double)(struc.testr[i])*as;
01890
01891 b += as;
01892 }
01893
01894 if( b <= 0. )
01895 {
01896 *tsqr = 0.;
01897 }
01898 else
01899 {
01900
01901 tmean = a/b;
01902 a = 0.;
01903
01904 ASSERT( nzone < struc.nzlim );
01905 for( i=0; i < nzone; i++ )
01906 {
01907 as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
01908 struc.ednstr[i];
01909 a += (POW2((double)(struc.testr[i]-tmean)))*as;
01910 }
01911 *tsqr = (realnum)(a/(b*(POW2(tmean))));
01912 }
01913
01914 return;
01915 }
01916
01917
01918 STATIC void gett2o3(realnum *tsqr)
01919 {
01920 long int i;
01921 double tmean;
01922 double a,
01923 as,
01924 b;
01925
01926 DEBUG_ENTRY( "gett2o3()" );
01927
01928 a = 0.;
01929 b = 0.;
01930 ASSERT(nzone<struc.nzlim);
01931 for( i=0; i < nzone; i++ )
01932 {
01933 as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
01934 (double)(struc.ednstr[i]);
01935 a += (double)(struc.testr[i])*as;
01936
01937
01938 b += as;
01939 }
01940
01941 if( b <= 0. )
01942 {
01943 *tsqr = 0.;
01944 }
01945
01946 else
01947 {
01948
01949 tmean = a/b;
01950 a = 0.;
01951 ASSERT( nzone < struc.nzlim );
01952 for( i=0; i < nzone; i++ )
01953 {
01954 as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
01955 struc.ednstr[i];
01956 a += (POW2((double)(struc.testr[i]-tmean)))*as;
01957 }
01958 *tsqr = (realnum)(a/(b*(POW2(tmean))));
01959 }
01960 return;
01961 }