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