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