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