cloudy trunk

prt_final.cpp

Go to the documentation of this file.
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 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated for cloudy by doxygen 1.7.3