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 Feb. 9, 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 "(Dere, K. P., Landi, E., Young, P. R., Del Zanna, G., Landini, M. & Mason, H. E. 2009, A&A, 498,\n 915) using version %s.\n",atmdat.chVersion);
00069 }
00070
00071 if( h2.lgH2ON )
00072 {
00073 fprintf( ioQQQ, " Much of the H_2 data is from "
00074 "Wrathmall & Flower 2007, JPhysB, 40 3221\n"
00075 " Abgrall et al, 1994, Can. J. Phys., 72, 856 at http://molat.obspm.fr"
00076 "/index.php?page=pages/Molecules/H2/H2can94.php\n");
00077 }
00078
00079 fprintf( ioQQQ, " Please cite those papers that played a significant role in this calculation.\n");
00080 }
00081
00082
00083 void PrintCenterLine(FILE* io,
00084 const char chLine[],
00085 size_t ArrLen,
00086 size_t LineLen)
00087 {
00088 unsigned long StrLen = min(strlen(chLine),ArrLen);
00089 ASSERT( StrLen < LineLen );
00090 unsigned long pad = (LineLen-StrLen)/2;
00091 for( unsigned long i=0; i < pad; ++i )
00092 fprintf( io, " " );
00093 fprintf( io, "%s\n", chLine );
00094 }
00095
00096
00097 STATIC void gett2o3(realnum *tsqr);
00098
00099
00100 STATIC void gett2(realnum *tsqr);
00101
00102
00103 inline void PrintRatio(double q1, double q2)
00104 {
00105 double ratio = ( q2 > SMALLFLOAT ) ? q1/q2 : 0.;
00106 fprintf( ioQQQ, " " );
00107 fprintf( ioQQQ, PrintEfmt("%9.2e", ratio) );
00108 return;
00109 }
00110
00111
00112 void PrtFinal(void)
00113 {
00114 short int *lgPrt;
00115 realnum *wavelength;
00116 realnum *sclsav , *scaled;
00117 long int *ipSortLines;
00118 double *xLog_line_lumin;
00119 char **chSLab;
00120 char *chSTyp;
00121 char chCKey[5],
00122 chGeo[7],
00123 chPlaw[21];
00124
00125 long int
00126 i,
00127 ipEmType ,
00128 ipNegIntensity[33],
00129 ip2500,
00130 ip2kev,
00131 iprnt,
00132 j,
00133 nline,
00134 nNegIntenLines;
00135 double o4363,
00136 bacobs,
00137 absint,
00138 bacthn,
00139 hbcab,
00140 hbeta,
00141 o5007;
00142
00143 double a,
00144 ajmass,
00145 ajmin,
00146 alfox,
00147 bot,
00148 bottom,
00149 bremtk,
00150 coleff,
00151
00152 d[8],
00153 dmean,
00154 ferode,
00155 he4686,
00156 he5876,
00157 heabun,
00158 hescal,
00159 pion,
00160 plow,
00161 powerl,
00162 qion,
00163 qlow,
00164 rate,
00165 ratio,
00166 reclin,
00167 rjeans,
00168 snorm,
00169 tmean,
00170 top,
00171 THI,
00172 uhel,
00173 uhl,
00174 usp,
00175 wmean,
00176 znit;
00177
00178 double bac,
00179 tel,
00180 x;
00181
00182 DEBUG_ENTRY( "PrtFinal()" );
00183
00184
00185
00186 if( !called.lgTalk || (lgAbort && nzone< 1) )
00187 {
00188 return;
00189 }
00190
00191
00192
00193
00194 ASSERT( LineSave.nsum > 1 );
00195
00196
00197 fprintf( ioQQQ, "\f\n");
00198 fprintf( ioQQQ, "%23c", ' ' );
00199 int len = (int)strlen(t_version::Inst().chVersion);
00200 int repeat = (72-len)/2;
00201 for( i=0; i < repeat; ++i )
00202 fprintf( ioQQQ, "*" );
00203 fprintf( ioQQQ, "> Cloudy %s <", t_version::Inst().chVersion );
00204 for( i=0; i < 72-repeat-len; ++i )
00205 fprintf( ioQQQ, "*" );
00206 fprintf( ioQQQ, "\n" );
00207
00208 for( i=0; i <= input.nSave; i++ )
00209 {
00210 char chCard[INPUT_LINE_LENGTH];
00211
00212
00213
00214 cap4(chCKey,input.chCardSav[i]);
00215
00216
00217 strcpy( chCard , input.chCardSav[i] );
00218 caps( chCard );
00219
00220
00221
00222 if( (strcmp(chCKey,"CONT")!= 0) && !nMatch( "HIDE" , chCard) )
00223 fprintf( ioQQQ, "%23c* %-80s*\n", ' ', input.chCardSav[i] );
00224 }
00225
00226
00227 if( rfield.uh > 0. )
00228 {
00229 a = log10(rfield.uh);
00230 }
00231 else
00232 {
00233 a = -37.;
00234 }
00235
00236 fprintf( ioQQQ,
00237 " *********************************> Log(U):%6.2f <*********************************\n",
00238 a );
00239
00240 if( t_version::Inst().nBetaVer > 0 )
00241 {
00242 fprintf( ioQQQ,
00243 "\n This is a beta test version of the code, and is intended for testing only.\n\n" );
00244 }
00245
00246 if( warnings.lgWarngs )
00247 {
00248 fprintf( ioQQQ, " \n" );
00249 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00250 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00251 fprintf( ioQQQ, " >>>>>>>>>> Warning! Warnings exist, this calculation has serious problems.\n" );
00252 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00253 fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
00254 fprintf( ioQQQ, " \n" );
00255 }
00256 else if( warnings.lgCautns )
00257 {
00258 fprintf( ioQQQ,
00259 " >>>>>>>>>> Cautions are present.\n" );
00260 }
00261
00262 if( conv.lgBadStop )
00263 {
00264 fprintf( ioQQQ,
00265 " >>>>>>>>>> The calculation stopped for unintended reasons.\n" );
00266 }
00267
00268 if( iterations.lgIterAgain )
00269 {
00270 fprintf( ioQQQ,
00271 " >>>>>>>>>> Another iteration is needed.\n" );
00272 }
00273
00274
00275 if( geometry.lgSphere )
00276 {
00277 strcpy( chGeo, "Closed" );
00278 }
00279 else
00280 {
00281 strcpy( chGeo, " Open" );
00282 }
00283
00284
00285 if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
00286 {
00287 strcpy( chPlaw, " Constant Pressure " );
00288 }
00289
00290 else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
00291 {
00292 strcpy( chPlaw, " Constant Density " );
00293 }
00294
00295 else if( (strcmp(dense.chDenseLaw,"POWD") == 0 || strcmp(dense.chDenseLaw
00296 ,"POWR") == 0) || strcmp(dense.chDenseLaw,"POWC") == 0 )
00297 {
00298 strcpy( chPlaw, " Power Law Density " );
00299 }
00300
00301 else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
00302 {
00303 strcpy( chPlaw, " Rapid Fluctuations " );
00304 }
00305
00306 else if( strncmp(dense.chDenseLaw , "DLW" , 3) == 0 )
00307 {
00308 strcpy( chPlaw, " Special Density Lw " );
00309 }
00310
00311 else if( strcmp(dense.chDenseLaw,"HYDR") == 0 )
00312 {
00313 strcpy( chPlaw, " Hydrostatic Equlib " );
00314 }
00315
00316 else if( strcmp(dense.chDenseLaw,"WIND") == 0 )
00317 {
00318 strcpy( chPlaw, " Radia Driven Wind " );
00319 }
00320
00321 else if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00322 {
00323 strcpy( chPlaw, " Globule " );
00324 }
00325
00326 else
00327 {
00328 strcpy( chPlaw, " UNRECOGNIZED CPRES " );
00329 }
00330
00331 fprintf( ioQQQ,
00332 "\n Emission Line Spectrum. %20.20sModel. %6.6s geometry. Iteration %ld of %ld.\n",
00333 chPlaw, chGeo, iteration, iterations.itermx + 1 );
00334
00335
00336 if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth )
00337 {
00338 char chLine[INPUT_LINE_LENGTH];
00339 const char* chUnit;
00340 if( geometry.iEmissPower == 1 && !geometry.lgSizeSet )
00341 chUnit = "/arcsec";
00342 else if( geometry.iEmissPower == 0 && !geometry.lgSizeSet )
00343 chUnit = "/arcsec^2";
00344 else
00345 chUnit = "";
00346
00347 sprintf( chLine, "Flux observed at the Earth (erg/s/cm^2%s).", chUnit );
00348 PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
00349 }
00350 else if( prt.lgSurfaceBrightness )
00351 {
00352 char chLine[INPUT_LINE_LENGTH];
00353 const char* chUnit;
00354 if( prt.lgSurfaceBrightness_SR )
00355 chUnit = "sr";
00356 else
00357 chUnit = "arcsec^2";
00358
00359 sprintf( chLine, "Surface brightness (erg/s/cm^2/%s).", chUnit );
00360 PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
00361 }
00362 else if( radius.lgPredLumin )
00363 {
00364 char chLine[INPUT_LINE_LENGTH];
00365 const char* chUnit;
00366 if( geometry.iEmissPower == 2 )
00367 chUnit = "erg/s";
00368 else if( geometry.iEmissPower == 1 )
00369 chUnit = "erg/s/cm";
00370 else if( geometry.iEmissPower == 0 )
00371 chUnit = "erg/s/cm^2";
00372 else
00373 TotalInsanity();
00374
00375 char chCoverage[INPUT_LINE_LENGTH];
00376 if( fp_equal( geometry.covgeo, realnum(1.) ) )
00377 sprintf( chCoverage, "with full coverage" );
00378 else
00379 sprintf( chCoverage, "with a covering factor of %.1f%%", geometry.covgeo*100. );
00380
00381 if( radius.lgCylnOn )
00382 sprintf( chLine, "Luminosity (%s) emitted by a partial cylinder %s.", chUnit, chCoverage );
00383 else
00384 sprintf( chLine, "Luminosity (%s) emitted by a shell %s.", chUnit, chCoverage );
00385 PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
00386
00387 if( geometry.iEmissPower != 2 )
00388 {
00389 const char* chAper;
00390 if( geometry.iEmissPower == 1 )
00391 chAper = "long slit";
00392 else if( geometry.iEmissPower == 0 )
00393 chAper = "pencil beam";
00394 else
00395 TotalInsanity();
00396
00397 sprintf( chLine, "Observed through a %s with aperture covering factor %.1f%%.",
00398 chAper, geometry.covaper*100. );
00399 PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
00400 }
00401 }
00402 else
00403 {
00404 char chLine[INPUT_LINE_LENGTH];
00405 sprintf( chLine, "Intensity (erg/s/cm^2)." );
00406 PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
00407 }
00408
00409 fprintf( ioQQQ, "\n" );
00410
00411
00412
00413
00414
00415
00416 if( !atmdat.lgHCaseBOK[1][ipHYDROGEN] )
00417 {
00418 static const int NWLH = 21;
00419
00420 realnum wlh[NWLH] = {6563.e0f, 4861.e0f, 4340.e0f, 4102.e0f, 1.875e4f, 1.282e4f,
00421 1.094e4f, 1.005e4f, 2.625e4f, 2.166e4f, 1.945e4f, 7.458e4f,
00422 4.653e4f, 3.740e4f, 4.051e4f, 7.458e4f, 3.296e4f, 1216.e0f,
00423 1026.e0f, 972.5e0f, 949.7e0f };
00424
00425
00426 for( i=0; i < LineSave.nsum; i++ )
00427 {
00428
00429
00430
00431 if( (strcmp( LineSv[i].chALab,"Ca B" )==0) ||
00432 (strcmp( LineSv[i].chALab,"Ca A" )==0) )
00433 {
00434 realnum errorwave;
00435
00436
00437
00438 errorwave = WavlenErrorGet( LineSv[i].wavelength );
00439 for( j=0; j<NWLH; ++j )
00440 {
00441 if( fabs(LineSv[i].wavelength-wlh[j] ) <= errorwave )
00442 {
00443 LineSv[i].SumLine[0] = 0.;
00444 LineSv[i].SumLine[1] = 0.;
00445 break;
00446 }
00447 }
00448 }
00449 }
00450 }
00451
00452 if( !atmdat.lgHCaseBOK[1][ipHELIUM] )
00453 {
00454
00455 static const int NWLHE = 20;
00456 realnum wlhe[NWLHE] = {1640.e0f, 1215.e0f, 1085.e0f, 1025.e0f, 4686.e0f, 3203.e0f,
00457 2733.e0f, 2511.e0f, 1.012e4f, 6560.e0f, 5412.e0f, 4860.e0f,
00458 1.864e4f, 1.163e4f, 9345.e0f, 8237.e0f, 303.8e0f, 256.3e0f,
00459 243.0e0f, 237.3e0f};
00460 for( i=0; i < LineSave.nsum; i++ )
00461 {
00462 if( (strcmp( LineSv[i].chALab,"Ca B" )==0) ||
00463 (strcmp( LineSv[i].chALab,"Ca A" )==0) )
00464 {
00465 realnum errorwave;
00466
00467
00468
00469 errorwave = WavlenErrorGet( LineSv[i].wavelength );
00470 for( j=0; j<NWLHE; ++j )
00471 {
00472 if( fabs(LineSv[i].wavelength-wlhe[j] ) <= errorwave )
00473 {
00474 LineSv[i].SumLine[0] = 0.;
00475 LineSv[i].SumLine[1] = 0.;
00476 break;
00477 }
00478 }
00479 }
00480 }
00481 }
00482
00483
00484
00485
00486
00487
00488
00489 if( cdLine("TOTL",4861.36f,&hbeta,&absint)<=0 )
00490 {
00491 if( dense.lgElmtOn[ipHYDROGEN] )
00492 {
00493
00494 fprintf( ioQQQ, " PrtFinal could not find TOTL 4861 with cdLine.\n" );
00495 cdEXIT(EXIT_FAILURE);
00496 }
00497 }
00498
00499 if( dense.lgElmtOn[ipHELIUM] )
00500 {
00501
00502
00503 if( cdLine("He 1",5875.61f,&he5876,&absint)<=0 )
00504 {
00505
00506 if( iso.numLevels_local[ipHE_LIKE][ipHELIUM] >= 14 )
00507 {
00508
00509 fprintf( ioQQQ, " PrtFinal could not find He 1 5876 with cdLine.\n" );
00510 cdEXIT(EXIT_FAILURE);
00511 }
00512 }
00513
00514
00515
00516 if( cdLine("He 2",4686.01f,&he4686,&absint)<=0 )
00517 {
00518
00519 if( iso.numLevels_local[ipH_LIKE][ipHELIUM] > 5 )
00520 {
00521
00522
00523 fprintf( ioQQQ, " PrtFinal could not find He 2 4686 with cdLine.\n" );
00524 cdEXIT(EXIT_FAILURE);
00525 }
00526 }
00527 }
00528 else
00529 {
00530 he5876 = 0.;
00531 absint = 0.;
00532 he4686 = 0.;
00533 }
00534
00535 if( hbeta > 0. )
00536 {
00537 heabun = (he4686*0.078 + he5876*0.739)/hbeta;
00538 }
00539 else
00540 {
00541 heabun = 0.;
00542 }
00543
00544 if( dense.lgElmtOn[ipHELIUM] )
00545 {
00546 hescal = heabun/(dense.gas_phase[ipHELIUM]/dense.gas_phase[ipHYDROGEN]);
00547 }
00548 else
00549 {
00550 hescal = 0.;
00551 }
00552
00553
00554
00555 if( cdLine("O 3",5007.,&o5007,&absint)<=0 )
00556 {
00557 if( dense.lgElmtOn[ipOXYGEN] )
00558 {
00559
00560 fprintf( ioQQQ, " PrtFinal could not find O 3 5007 with cdLine.\n" );
00561 cdEXIT(EXIT_FAILURE);
00562 }
00563 }
00564
00565 if( cdLine("TOTL",4363.,&o4363,&absint)<=0 )
00566 {
00567 if( dense.lgElmtOn[ipOXYGEN] )
00568 {
00569
00570 fprintf( ioQQQ, " PrtFinal could not find TOTL 4363 with cdLine.\n" );
00571 cdEXIT(EXIT_FAILURE);
00572 }
00573 }
00574
00575
00576 if( (o4363 != 0.) && (o5007 != 0.) )
00577 {
00578
00579
00580 bot = o5007 - o4363/oxy.o3enro;
00581
00582 if( bot > 0. )
00583 {
00584 ratio = o4363/bot;
00585 }
00586 else
00587 {
00588 ratio = 0.178;
00589 }
00590
00591 if( ratio > 0.177 )
00592 {
00593
00594 peimbt.toiiir = 1e10;
00595 }
00596 else
00597 {
00598
00599
00600
00601
00602
00603 oxy.o3cs12 = 2.2347f;
00604 oxy.o3cs13 = 0.29811f;
00605 ratio = ratio/1.3333/(oxy.o3cs13/oxy.o3cs12);
00606
00607 ratio = ratio/oxy.o3enro/oxy.o3br32;
00608 peimbt.toiiir = (realnum)(-oxy.o3ex23/log(ratio));
00609 }
00610 }
00611
00612 else
00613 {
00614 peimbt.toiiir = 0.;
00615 }
00616
00617 if( geometry.iEmissPower == 2 )
00618 {
00619
00620 if( cdLine("Bac ",3646.,&bacobs,&absint)<=0 )
00621 {
00622 fprintf( ioQQQ, " PrtFinal could not find Bac 3546 with cdLine.\n" );
00623 cdEXIT(EXIT_FAILURE);
00624 }
00625
00626
00627 if( hbeta > 0. )
00628 {
00629 bac = bacobs/hbeta;
00630 }
00631 else
00632 {
00633 bac = 0.;
00634 }
00635 }
00636 else
00637 {
00638 bac = 0.;
00639 }
00640
00641 if( bac > 0. )
00642 {
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661 x = 1.e0/log10(bac);
00662 tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 +
00663 x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));
00664
00665 if( tel > 1. && tel < 5. )
00666 {
00667 peimbt.tbac = (realnum)pow(10.,tel);
00668 }
00669 else
00670 {
00671 peimbt.tbac = 0.;
00672 }
00673 }
00674 else
00675 {
00676 peimbt.tbac = 0.;
00677 }
00678
00679 if( geometry.iEmissPower == 2 )
00680 {
00681
00682 if( cdLine("thin",3646.,&bacthn,&absint)<=0 )
00683 {
00684 fprintf( ioQQQ, " PrtFinal could not find thin 3646 with cdLine.\n" );
00685 cdEXIT(EXIT_FAILURE);
00686 }
00687
00688
00689 if( cdLine("Ca B",4861.36f,&hbcab,&absint)<=0 )
00690 {
00691 fprintf( ioQQQ, " PrtFinal could not find Ca B 4861 with cdLine.\n" );
00692 cdEXIT(EXIT_FAILURE);
00693 }
00694
00695 if( hbcab > 0. )
00696 {
00697 bacthn /= hbcab;
00698 }
00699 else
00700 {
00701 bacthn = 0.;
00702 }
00703 }
00704 else
00705 {
00706 bacthn = 0.;
00707 }
00708
00709 if( bacthn > 0. )
00710 {
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729 x = 1.e0/log10(bacthn);
00730 tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 +
00731 x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));
00732
00733 if( tel > 1. && tel < 5. )
00734 {
00735 peimbt.tbcthn = (realnum)pow(10.,tel);
00736 }
00737 else
00738 {
00739 peimbt.tbcthn = 0.;
00740 }
00741 }
00742 else
00743 {
00744 peimbt.tbcthn = 0.;
00745 }
00746
00747
00748
00749
00750 peimbt.tohyox = (realnum)((8.5*peimbt.toiiir - 7.5*peimbt.tbcthn - 228200. +
00751 sqrt(POW2(8.5*peimbt.toiiir-7.5*peimbt.tbcthn-228200.)+9.128e5*
00752 peimbt.tbcthn))/2.);
00753
00754 if( peimbt.tohyox > 0. )
00755 {
00756 peimbt.t2hyox = (realnum)((peimbt.tohyox - peimbt.tbcthn)/(1.7*peimbt.tohyox));
00757 }
00758 else
00759 {
00760 peimbt.t2hyox = 0.;
00761 }
00762
00763
00764
00765
00766
00767
00768 scaled = (realnum *)MALLOC( sizeof(realnum)*(unsigned long)LineSave.nsum);
00769
00770
00771 xLog_line_lumin = (double *)MALLOC( sizeof(double)*(unsigned long)LineSave.nsum);
00772
00773
00774
00775 lgPrt = (short int *)MALLOC( sizeof(short int)*(unsigned long)LineSave.nsum);
00776
00777
00778 wavelength = (realnum *)MALLOC( sizeof(realnum)*(unsigned long)LineSave.nsum);
00779
00780
00781 sclsav = (realnum *)MALLOC( sizeof(realnum)*(unsigned long)LineSave.nsum );
00782
00783
00784 chSTyp = (char *)MALLOC( sizeof(char)*(unsigned long)LineSave.nsum );
00785
00786
00787
00788
00789 chSLab = ((char**)MALLOC((unsigned long)LineSave.nsum*sizeof(char*)));
00790
00791
00792 for( i=0; i<LineSave.nsum; ++i)
00793 {
00794 chSLab[i] = (char*)MALLOC(5*sizeof(char));
00795 }
00796
00797
00798 ipSortLines = (long *)MALLOC( sizeof(long)*(unsigned long)LineSave.nsum );
00799
00800 ASSERT( LineSave.ipNormWavL >= 0 );
00801
00802
00803
00804 int nEmType = 2;
00805 if( prt.lgPrintLineCumulative && iteration > dynamics.n_initial_relax )
00806 nEmType = 4;
00807
00808 for( ipEmType=0; ipEmType<nEmType; ++ipEmType )
00809 {
00810
00811 snorm = LineSv[LineSave.ipNormWavL].SumLine[ipEmType];
00812
00813
00814 if( ((snorm <= SMALLDOUBLE ) || (LineSave.ipNormWavL < 0)) || (LineSave.ipNormWavL > LineSave.nsum) )
00815 {
00816 fprintf( ioQQQ, "\n\n >>PROBLEM Normalization line has small or zero intensity, its value was %.2e and its intensity was set to 1."
00817 "\n >>Please consider using another normalization line (this is set with the NORMALIZE command).\n" , snorm);
00818 fprintf( ioQQQ, " >>The relative intensities will be meaningless, and many lines may appear too faint.\n" );
00819 snorm = 1.;
00820 }
00821 for( i=0; i < LineSave.nsum; i++ )
00822 {
00823
00824
00825
00826
00827
00828 double scale = LineSv[i].SumLine[ipEmType]/snorm*LineSave.ScaleNormLine;
00829
00830 scale = MIN2(BIGFLOAT , scale );
00831 scale = MAX2( -BIGFLOAT , scale );
00832
00833
00834 scaled[i] = (realnum)scale;
00835
00836 if( LineSv[i].SumLine[ipEmType] > 0. )
00837 {
00838 xLog_line_lumin[i] = log10(LineSv[i].SumLine[ipEmType]) + radius.Conv2PrtInten;
00839 }
00840 else
00841 {
00842 xLog_line_lumin[i] = -38.;
00843 }
00844 }
00845
00846
00847 for( i=0; i < LineSave.nsum; i++ )
00848 {
00849
00850 if( (strcmp(LineSv[i].chALab,"Unit") == 0) || (strcmp(LineSv[i].chALab,"UntD") == 0) )
00851 lgPrt[i] = false;
00852 else if( strcmp(LineSv[i].chALab,"Coll") == 0 && !prt.lgPrnColl )
00853 lgPrt[i] = false;
00854 else if( strcmp(LineSv[i].chALab,"Pump") == 0 && !prt.lgPrnPump )
00855 lgPrt[i] = false;
00856 else if( strncmp(LineSv[i].chALab,"Inw",3) == 0 && !prt.lgPrnInwd )
00857 lgPrt[i] = false;
00858 else if( strcmp(LineSv[i].chALab,"Heat") == 0 && !prt.lgPrnHeat )
00859 lgPrt[i] = false;
00860
00861
00862 else if( !prt.lgPrnDiff &&
00863 ( (strcmp(LineSv[i].chALab,"nFnu") == 0) || (strcmp(LineSv[i].chALab,"nInu") == 0) ) )
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.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 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.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-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_itrzn());
01450
01451
01452 fprintf(ioQQQ, " File Opacity: F" );
01453
01454
01455 {
01456
01457 double xmass = log10( SDIV(dense.xMassTotal) );
01458 xmass += (realnum)(1.0992099 + 2.*log10(radius.rinner));
01459 fprintf(ioQQQ," Mass tot %.3f",
01460 xmass);
01461 }
01462 fprintf(ioQQQ,"\n");
01463
01464
01465
01466
01467 if( cdTemp( "opti",0,&THI,"volume" ) )
01468 {
01469 fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm opt.\n");
01470 TotalInsanity();
01471 }
01472 fprintf( ioQQQ, " Temps(21 cm) T(21cm/Ly a) ");
01473 PrintE82(ioQQQ, THI );
01474
01475
01476
01477
01478
01479 if( cdTemp( "21cm",0,&THI,"radius" ) )
01480 {
01481 fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm temp.\n");
01482 TotalInsanity();
01483 }
01484 fprintf(ioQQQ, " T(<nH/Tkin>) ");
01485 PrintE82(ioQQQ, THI);
01486
01487
01488
01489 if( cdTemp( "spin",0,&THI,"radius" ) )
01490 {
01491 fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm spin.\n");
01492 TotalInsanity();
01493 }
01494 fprintf(ioQQQ, " T(<nH/Tspin>) ");
01495 PrintE82(ioQQQ, THI);
01496
01497
01498 THI *= (1. - sexp( HFLines[0].Emis->TauIn ) );
01499 fprintf( ioQQQ, " TB21cm:");
01500 PrintE82(ioQQQ, THI);
01501 fprintf( ioQQQ, "\n");
01502
01503 fprintf( ioQQQ, " N(H0/Tspin) ");
01504 PrintE82(ioQQQ, colden.H0_ov_Tspin );
01505 fprintf( ioQQQ, " N(OH/Tkin) ");
01506 PrintE82(ioQQQ, colden.OH_ov_Tspin );
01507
01508
01509 bot = cdB21cm();
01510 fprintf( ioQQQ, " B(21cm) ");
01511 PrintE82(ioQQQ, bot );
01512
01513
01514 fprintf(ioQQQ, "\n");
01515
01516
01517 rate = timesc.TimeErode*2e-26;
01518 if( rate > SMALLFLOAT )
01519 {
01520 ferode = 1./rate;
01521 }
01522 else
01523 {
01524 ferode = 0.;
01525 }
01526
01527
01528 if( wind.acldr > 0. )
01529 {
01530 wind.AccelAver /= wind.acldr;
01531 }
01532 else
01533 {
01534 wind.AccelAver = 0.;
01535 }
01536
01537
01538 wmean = colden.wmas/SDIV(colden.TotMassColl);
01539
01540 dmean = colden.TotMassColl/SDIV(radius.depth_x_fillfac);
01541 tmean = colden.tmas/SDIV(colden.TotMassColl);
01542
01543 wmean = colden.wmas/SDIV(colden.TotMassColl);
01544
01545 fprintf( ioQQQ, " <a>:");
01546 PrintE82(ioQQQ , wind.AccelAver);
01547 fprintf( ioQQQ, " erdeFe");
01548 PrintE71(ioQQQ , ferode);
01549 fprintf( ioQQQ, " Tcompt");
01550 PrintE82(ioQQQ , timesc.tcmptn);
01551 fprintf( ioQQQ, " Tthr");
01552 PrintE82(ioQQQ , timesc.time_therm_long);
01553 fprintf( ioQQQ, " <Tden>: ");
01554 PrintE82(ioQQQ , tmean);
01555 fprintf( ioQQQ, " <dens>:");
01556 PrintE82(ioQQQ , dmean);
01557 fprintf( ioQQQ, " <MolWgt>");
01558 PrintE82(ioQQQ , wmean);
01559 fprintf(ioQQQ,"\n");
01560
01561
01562 if( tmean > 0. )
01563 {
01564 rjeans = 7.79637 + (log10(tmean) - log10(dense.wmole) - log10(dense.xMassDensity*
01565 geometry.FillFac))/2.;
01566 }
01567 else
01568 {
01569
01570 rjeans = 40.f;
01571 }
01572
01573 if( rjeans < 36. )
01574 {
01575 rjeans = (double)pow(10.,rjeans);
01576
01577 ajmass = 3.*log10(rjeans/2.) + log10(4.*PI/3.*dense.xMassDensity*
01578 geometry.FillFac) - log10(SOLAR_MASS);
01579 if( ajmass < 37 )
01580 {
01581 ajmass = pow(10.,ajmass);
01582 }
01583 else
01584 {
01585 ajmass = 0.;
01586 }
01587 }
01588 else
01589 {
01590 rjeans = 0.;
01591 ajmass = 0.;
01592 }
01593
01594
01595 ajmin = colden.ajmmin - log10(SOLAR_MASS);
01596 if( ajmin < 37 )
01597 {
01598 ajmin = pow(10.,ajmin);
01599 }
01600 else
01601 {
01602 ajmin = 0.;
01603 }
01604
01605
01606 if( rfield.anu[rfield.nflux-1] > 150. )
01607 {
01608
01609 ip2kev = ipoint(147.);
01610 ip2500 = ipoint(0.3648);
01611
01612
01613 bottom = rfield.flux[0][ip2500-1]*
01614 rfield.anu[ip2500-1]/rfield.widflx[ip2500-1];
01615
01616 top = rfield.flux[0][ip2kev-1]*
01617 rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1];
01618
01619
01620 if( bottom > 1e-30 && top > 1e-30 )
01621 {
01622 ratio = log10(top) - log10(bottom);
01623 if( ratio < 36. && ratio > -36. )
01624 {
01625 ratio = pow(10.,ratio);
01626 }
01627 else
01628 {
01629 ratio = 0.;
01630 }
01631 }
01632
01633 else
01634 {
01635 ratio = 0.;
01636 }
01637
01638 if( ratio > 0. )
01639 {
01640
01641 double freq_ratio = rfield.anu[ip2kev-1]/rfield.anu[ip2500-1];
01642 alfox = log(ratio)/log(freq_ratio);
01643 }
01644 else
01645 {
01646 alfox = 0.;
01647 }
01648 }
01649 else
01650 {
01651 alfox = 0.;
01652 }
01653
01654 if( colden.rjnmin < 37 )
01655 {
01656 colden.rjnmin = (realnum)pow((realnum)10.f,colden.rjnmin);
01657 }
01658 else
01659 {
01660 colden.rjnmin = FLT_MAX;
01661 }
01662
01663 fprintf( ioQQQ, " Mean Jeans l(cm)");
01664 PrintE82(ioQQQ, rjeans);
01665 fprintf( ioQQQ, " M(sun)");
01666 PrintE82(ioQQQ, ajmass);
01667 fprintf( ioQQQ, " smallest: len(cm):");
01668 PrintE82(ioQQQ, colden.rjnmin);
01669 fprintf( ioQQQ, " M(sun):");
01670 PrintE82(ioQQQ, ajmin);
01671 fprintf( ioQQQ, " a_ox tran:%6.2f\n" , alfox);
01672
01673 fprintf( ioQQQ, " Rion:");
01674 double R_ion;
01675 if( rfield.lgUSphON )
01676 R_ion = rfield.rstrom;
01677 else
01678 R_ion = radius.Radius;
01679 PrintE93(ioQQQ, R_ion);
01680 fprintf( ioQQQ, " Dist:");
01681 PrintE82(ioQQQ, radius.distance);
01682 fprintf( ioQQQ, " Diam:");
01683 if( radius.distance > 0. )
01684 PrintE93(ioQQQ, 2.*R_ion*AS1RAD/radius.distance);
01685 else
01686 PrintE93(ioQQQ, 0.);
01687 fprintf( ioQQQ, "\n");
01688
01689 if( prt.lgPrintTime )
01690 {
01691
01692 alfox = cdExecTime();
01693 }
01694 else
01695 {
01696
01697
01698 alfox = 0.;
01699 }
01700
01701
01702 fprintf( ioQQQ,
01703 " Hatom level%3ld HatomType:%4.4s HInducImp %2c"
01704 " He tot level:%3ld He2 level: %3ld ExecTime%8.1f\n",
01705
01706 iso.numLevels_local[ipH_LIKE][ipHYDROGEN],
01707 hydro.chHTopType,
01708 TorF(hydro.lgHInducImp),
01709
01710 iso.numLevels_local[ipHE_LIKE][ipHELIUM],
01711
01712 iso.numLevels_local[ipH_LIKE][ipHELIUM] ,
01713 alfox );
01714
01715
01716 fprintf( ioQQQ,
01717 " ConvrgError(%%) <eden>%7.3f MaxEden%7.3f <H-C>%7.2f Max(H-C)%8.2f <Press>%8.3f MaxPrs er%7.3f\n",
01718 conv.AverEdenError/nzone*100. ,
01719 conv.BigEdenError*100. ,
01720 conv.AverHeatCoolError/nzone*100. ,
01721 conv.BigHeatCoolError*100. ,
01722 conv.AverPressError/nzone*100. ,
01723 conv.BigPressError*100. );
01724
01725 fprintf(ioQQQ,
01726 " Continuity(%%) chng Te%6.1f elec den%6.1f n(H2)%7.1f n(CO) %7.1f\n",
01727 phycon.BigJumpTe*100. ,
01728 phycon.BigJumpne*100. ,
01729 phycon.BigJumpH2*100. ,
01730 phycon.BigJumpCO*100. );
01731
01732
01733 fprintf( ioQQQ, "\n Averaged Quantities\n" );
01734 fprintf( ioQQQ, " Te Te(Ne) Te(NeNp) Te(NeHe+)Te(NeHe2+) Te(NeO+) Te(NeO2+)"
01735 " Te(H2) N(H) Ne(O2+) Ne(Np)\n" );
01736 static const char* weight[3] = { "Radius", "Area", "Volume" };
01737 int dmax = geometry.lgGeoPP ? 1 : 3;
01738 for( int dim=0; dim < dmax; ++dim )
01739 {
01740 fprintf( ioQQQ, " %6s:", weight[dim] );
01741
01742 PrintRatio( mean.TempMean[dim][0], mean.TempMean[dim][1] );
01743
01744 PrintRatio( mean.TempEdenMean[dim][0], mean.TempEdenMean[dim][1] );
01745
01746 PrintRatio( mean.TempIonEdenMean[dim][ipHYDROGEN][1][0], mean.TempIonEdenMean[dim][ipHYDROGEN][1][1] );
01747 PrintRatio( mean.TempIonEdenMean[dim][ipHELIUM][1][0], mean.TempIonEdenMean[dim][ipHELIUM][1][1] );
01748 PrintRatio( mean.TempIonEdenMean[dim][ipHELIUM][2][0], mean.TempIonEdenMean[dim][ipHELIUM][2][1] );
01749 PrintRatio( mean.TempIonEdenMean[dim][ipOXYGEN][1][0], mean.TempIonEdenMean[dim][ipOXYGEN][1][1] );
01750 PrintRatio( mean.TempIonEdenMean[dim][ipOXYGEN][2][0], mean.TempIonEdenMean[dim][ipOXYGEN][2][1] );
01751
01752 PrintRatio( mean.TempIonMean[dim][ipHYDROGEN][2][0], mean.TempIonMean[dim][ipHYDROGEN][2][1] );
01753
01754 PrintRatio( mean.xIonMean[dim][ipHYDROGEN][0][1], mean.TempMean[dim][1] );
01755
01756 PrintRatio( mean.TempIonEdenMean[dim][ipOXYGEN][2][1], mean.TempIonMean[dim][ipOXYGEN][2][1] );
01757 PrintRatio( mean.TempIonEdenMean[dim][ipHYDROGEN][1][1], mean.TempIonMean[dim][ipHYDROGEN][1][1] );
01758 fprintf( ioQQQ, "\n" );
01759 }
01760
01761
01762
01763 if( dense.gas_phase[ipHYDROGEN] < peimbt.tsqden )
01764 {
01765 fprintf( ioQQQ, " \n" );
01766
01767
01768 fprintf( ioQQQ, " Peimbert T(OIIIr)");
01769 PrintE82( ioQQQ , peimbt.toiiir);
01770
01771
01772 fprintf( ioQQQ, " T(Bac)");
01773 PrintE82( ioQQQ , peimbt.tbac);
01774
01775
01776 fprintf( ioQQQ, " T(Hth)");
01777 PrintE82( ioQQQ , peimbt.tbcthn);
01778
01779
01780 fprintf( ioQQQ, " t2(Hstrc)");
01781 fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hstr));
01782
01783
01784 fprintf( ioQQQ, " T(O3-BAC)");
01785 PrintE82( ioQQQ , peimbt.tohyox);
01786
01787
01788 fprintf( ioQQQ, " t2(O3-BC)");
01789 fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hyox));
01790
01791
01792 fprintf( ioQQQ, " t2(O3str)");
01793 fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2o3str));
01794
01795 fprintf( ioQQQ, "\n");
01796
01797 if( gv.lgDustOn() )
01798 {
01799 fprintf( ioQQQ, " Be careful: grains exist. This spectrum was not corrected for reddening before analysis.\n" );
01800 }
01801 }
01802
01803
01804
01805 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
01806 {
01807 char chQHeat;
01808 double AV , AB;
01809 double total_dust2gas = 0.;
01810
01811 fprintf( ioQQQ, "\n Average Grain Properties (over radius):\n" );
01812
01813 for( size_t i0=0; i0 < gv.bin.size(); i0 += 10 )
01814 {
01815 if( i0 > 0 )
01816 fprintf( ioQQQ, "\n" );
01817
01818
01819 size_t i1 = min(i0+10,gv.bin.size());
01820
01821 fprintf( ioQQQ, " " );
01822 for( size_t nd=i0; nd < i1; nd++ )
01823 {
01824 chQHeat = (char)(gv.bin[nd]->lgEverQHeat ? '*' : ' ');
01825 fprintf( ioQQQ, "%-12.12s%c", gv.bin[nd]->chDstLab, chQHeat );
01826 }
01827 fprintf( ioQQQ, "\n" );
01828
01829 fprintf( ioQQQ, " nd:" );
01830 for( size_t nd=i0; nd < i1; nd++ )
01831 {
01832 if( nd != i0 ) fprintf( ioQQQ," " );
01833 fprintf( ioQQQ, "%7ld ", (unsigned long)nd );
01834 }
01835 fprintf( ioQQQ, "\n" );
01836
01837 fprintf( ioQQQ, " <Tgr>:" );
01838 for( size_t nd=i0; nd < i1; nd++ )
01839 {
01840 if( nd != i0 ) fprintf( ioQQQ," " );
01841 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdust/radius.depth_x_fillfac));
01842 }
01843 fprintf( ioQQQ, "\n" );
01844
01845 fprintf( ioQQQ, " <Vel>:" );
01846 for( size_t nd=i0; nd < i1; nd++ )
01847 {
01848 if( nd != i0 ) fprintf( ioQQQ," " );
01849 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdft/radius.depth_x_fillfac));
01850 }
01851 fprintf( ioQQQ, "\n" );
01852
01853 fprintf( ioQQQ, " <Pot>:" );
01854 for( size_t nd=i0; nd < i1; nd++ )
01855 {
01856 if( nd != i0 ) fprintf( ioQQQ," " );
01857 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdpot/radius.depth_x_fillfac));
01858 }
01859 fprintf( ioQQQ, "\n" );
01860
01861 fprintf( ioQQQ, " <D/G>:" );
01862 for( size_t nd=i0; nd < i1; nd++ )
01863 {
01864 if( nd != i0 ) fprintf( ioQQQ," " );
01865 fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avDGRatio/radius.depth_x_fillfac));
01866
01867 total_dust2gas += gv.bin[nd]->avDGRatio/radius.depth_x_fillfac;
01868 }
01869 fprintf( ioQQQ, "\n" );
01870 }
01871
01872 fprintf(ioQQQ," Dust to gas ratio (by mass):");
01873 fprintf(ioQQQ,PrintEfmt("%10.3e", total_dust2gas));
01874
01875
01876
01877
01878
01879 AV = rfield.extin_mag_V_point/SDIV(colden.colden[ipCOL_HTOT]);
01880 AB = rfield.extin_mag_B_point/SDIV(colden.colden[ipCOL_HTOT]);
01881
01882 fprintf(ioQQQ,", A(V)/N(H)(pnt):%.3e, (ext):%.3e",
01883 AV,
01884 rfield.extin_mag_V_extended/SDIV(colden.colden[ipCOL_HTOT]) );
01885
01886
01887 fprintf(ioQQQ,", R:");
01888
01889
01890 if( fabs(AB-AV)>SMALLFLOAT )
01891 {
01892 fprintf(ioQQQ,"%.3e", AV/(AB-AV) );
01893 }
01894 else
01895 {
01896 fprintf(ioQQQ,"%.3e", 0. );
01897 }
01898 fprintf(ioQQQ," AV(ext):%.3e (pnt):%.3e\n",
01899 rfield.extin_mag_V_extended,
01900 rfield.extin_mag_V_point);
01901 }
01902
01903
01904 free( wavelength );
01905 free( ipSortLines );
01906 free( sclsav );
01907 free( lgPrt );
01908 free( chSTyp );
01909
01910
01911 for( i=0; i < LineSave.nsum; ++i )
01912 {
01913 free( chSLab[i] );
01914 }
01915 free( chSLab );
01916
01917 free( scaled );
01918 free( xLog_line_lumin );
01919
01920
01921 if( !prt.lgPrtShort && called.lgTalk )
01922 {
01923
01924
01925 PrtAllTau();
01926
01927
01928 if( iterations.lgLastIt )
01929 {
01930
01931 if( prt.lgPrintColumns )
01932 PrtColumns(ioQQQ);
01933
01934 PrtMeanIon('i', false, ioQQQ);
01935
01936 PrtMeanIon('i', true , ioQQQ);
01937
01938 PrtMeanIon('t', false , ioQQQ);
01939
01940 PrtMeanIon('t', true , ioQQQ);
01941
01942 PrtContinuum();
01943 }
01944 }
01945
01946
01947 fprintf( ioQQQ, "%s\n\n", input.chTitle );
01948 return;
01949 }
01950
01951
01952
01953 long int StuffComment( const char * chComment )
01954 {
01955 long int n , i;
01956
01957 DEBUG_ENTRY( "StuffComment()" );
01958
01959
01960 if( LineSave.ipass == 0 )
01961 {
01962 if( LineSave.nComment >= NHOLDCOMMENTS )
01963 {
01964 fprintf( ioQQQ, " Too many comments have been entered into line array; increase the value of NHOLDCOMMENTS.\n" );
01965 cdEXIT(EXIT_FAILURE);
01966 }
01967
01968
01969 static const int NWIDTH = 33;
01970 strcpy( LineSave.chHoldComments[LineSave.nComment], chComment );
01971
01972
01973 n = (long)strlen( LineSave.chHoldComments[LineSave.nComment] );
01974 for( i=0; i<NWIDTH-n-1-6; ++i )
01975 {
01976 strcat( LineSave.chHoldComments[LineSave.nComment], ".");
01977 }
01978
01979 strcat( LineSave.chHoldComments[LineSave.nComment], "..");
01980
01981 for( i=0; i<6; ++i )
01982 {
01983 strcat( LineSave.chHoldComments[LineSave.nComment], " ");
01984 }
01985 }
01986
01987 ++LineSave.nComment;
01988 return( LineSave.nComment-1 );
01989 }
01990
01991
01992 STATIC void gett2(realnum *tsqr)
01993 {
01994 long int i;
01995
01996 double tmean;
01997 double a,
01998 as,
01999 b;
02000
02001 DEBUG_ENTRY( "gett2()" );
02002
02003
02004 a = 0.;
02005 b = 0.;
02006
02007 ASSERT( nzone < struc.nzlim);
02008
02009
02010
02011 for( i=0; i < nzone; i++ )
02012 {
02013 as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
02014 (double)(struc.ednstr[i]);
02015 a += (double)(struc.testr[i])*as;
02016
02017 b += as;
02018 }
02019
02020 if( b <= 0. )
02021 {
02022 *tsqr = 0.;
02023 }
02024 else
02025 {
02026
02027 tmean = a/b;
02028 a = 0.;
02029
02030 ASSERT( nzone < struc.nzlim );
02031 for( i=0; i < nzone; i++ )
02032 {
02033 as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
02034 struc.ednstr[i];
02035 a += (POW2((double)(struc.testr[i]-tmean)))*as;
02036 }
02037 *tsqr = (realnum)(a/(b*(POW2(tmean))));
02038 }
02039
02040 return;
02041 }
02042
02043
02044 STATIC void gett2o3(realnum *tsqr)
02045 {
02046 long int i;
02047 double tmean;
02048 double a,
02049 as,
02050 b;
02051
02052 DEBUG_ENTRY( "gett2o3()" );
02053
02054 a = 0.;
02055 b = 0.;
02056 ASSERT(nzone<struc.nzlim);
02057
02058
02059
02060 for( i=0; i < nzone; i++ )
02061 {
02062 as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
02063 (double)(struc.ednstr[i]);
02064 a += (double)(struc.testr[i])*as;
02065
02066
02067 b += as;
02068 }
02069
02070 if( b <= 0. )
02071 {
02072 *tsqr = 0.;
02073 }
02074
02075 else
02076 {
02077
02078 tmean = a/b;
02079 a = 0.;
02080 ASSERT( nzone < struc.nzlim );
02081 for( i=0; i < nzone; i++ )
02082 {
02083 as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
02084 struc.ednstr[i];
02085 a += (POW2((double)(struc.testr[i]-tmean)))*as;
02086 }
02087 *tsqr = (realnum)(a/(b*(POW2(tmean))));
02088 }
02089 return;
02090 }