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 /*CoolSum total cooling from all entries into cooling stack */ 00004 /*CoolZero set cooling and heating stack to zero */ 00005 /*CoolAdd add coolants to the cooling stack, called in evaluation of cooling function */ 00006 #include "cddefines.h" 00007 #include "taulines.h" 00008 #include "lines_service.h" 00009 #include "thermal.h" 00010 #include "cooling.h" 00011 00012 /*CoolAdd add coolants to the cooling stack, called in evaluation of cooling function */ 00013 void CoolAdd( 00014 const char *chLabel, 00015 realnum lambda, 00016 double cool) 00017 { 00018 00019 DEBUG_ENTRY( "CoolAdd()" ); 00020 00021 /* this flag indicates (true) that we are between when cooling was set to 00022 * zero with call to CoolZero, and when final sum was used. Any call 00023 * after final summation (false) will be ignored and so is fatal error */ 00024 ASSERT( thermal.lgCoolEvalOK ); 00025 00026 /* this can be done with an assert since these results cannot possibly 00027 * depend on user input */ 00028 ASSERT( thermal.ncltot < NCOLNT ); 00029 00030 /* copy coolant label into stack */ 00031 ASSERT( strlen( thermal.chClntLab[thermal.ncltot]) < NCOLNT_LAB_LEN ); 00032 strcpy( thermal.chClntLab[thermal.ncltot], chLabel); 00033 00034 /* now the wavelength */ 00035 thermal.collam[thermal.ncltot] = lambda; 00036 00037 /* normal line cooling */ 00038 thermal.cooling[thermal.ncltot] = MAX2(0.,cool); 00039 00040 /* possible line heating - not supposed to be done this way! 00041 * this is intrinsic positive number, to be added to heating */ 00042 thermal.heatnt[thermal.ncltot] = MAX2(0.,-cool); 00043 00044 /* now increment counter, this is the number of coolants entered */ 00045 thermal.ncltot += 1; 00046 return; 00047 } 00048 00049 /*CoolZero set cooling and heating stack to zero */ 00050 void CoolZero(void) 00051 { 00052 00053 DEBUG_ENTRY( "CoolZero()" ); 00054 00055 thermal.ncltot = 0; 00056 thermal.dCooldT = 0.; 00057 00058 /* >>chng 03 nov 29, from explicit loop to memset to save time */ 00059 memset(thermal.cooling , 0 , NCOLNT*sizeof(thermal.cooling[0] ) ); 00060 memset(thermal.heatnt , 0 , NCOLNT*sizeof(thermal.heatnt[0] ) ); 00061 00062 /* this flag indicates that it is ok to add coolants to cooling 00063 * stack since between first zero, and final sum - CoolAdd checks 00064 * that this is true */ 00065 thermal.lgCoolEvalOK = true; 00066 return; 00067 } 00068 00069 /*CoolSum total cooling from all entries into cooling stack */ 00070 void CoolSum(double *total) 00071 { 00072 long int i; 00073 00074 DEBUG_ENTRY( "CoolSum()" ); 00075 00076 /* routine to add together all line heating and cooling */ 00077 00078 *total = 0.; 00079 thermal.coolheat = 0.; 00080 /* this is sum of agents that should be coolants 00081 * coolheat will be coolants that came out as heat */ 00082 for( i=0; i < thermal.ncltot; i++ ) 00083 { 00084 *total += thermal.cooling[i]; 00085 thermal.coolheat += thermal.heatnt[i]; 00086 } 00087 thermal.heating[0][12] = thermal.coolheat; 00088 00089 /* make comment if negative cooling ever significant */ 00090 if( thermal.htot > 0. ) 00091 { 00092 if( thermal.coolheat/thermal.htot > 0.01 ) 00093 { 00094 /* CoolHeatMax was set to zero at start of calc, we want very biggest */ 00095 for( i=0; i < thermal.ncltot; i++ ) 00096 { 00097 if( thermal.heatnt[i]/thermal.htot > thermal.CoolHeatMax ) 00098 { 00099 thermal.CoolHeatMax = (realnum)(thermal.heatnt[i]/thermal.htot); 00100 thermal.wlCoolHeatMax = thermal.collam[i]; 00101 strcpy( thermal.chCoolHeatMax, thermal.chClntLab[i] ); 00102 } 00103 } 00104 } 00105 } 00106 00107 /* this sum of lines that were heat sources - this 00108 * part was not counted as heating in call to cooling add routine 00109 * since atom_level2 and atom_level3 cooling routines separate this out 00110 * into ->cool and ->heat - this does 00111 * NOT double count line heating */ 00112 00113 thermal.heatl = 0.; 00114 for( i=0; i < nWindLine; i++ ) 00115 { 00116 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00117 thermal.heatl += TauLine2[i].Coll.heat; 00118 } 00119 00120 for( i=1; i <= nLevel1; i++ ) 00121 { 00122 thermal.heatl += TauLines[i].Coll.heat; 00123 } 00124 00125 /* the Chianti and Leiden lines */ 00126 for( i=0; i < linesAdded2; i++ ) 00127 { 00128 thermal.heatl += atmolEmis[i].tran->Coll.heat; 00129 } 00130 00131 # if 0 00132 /* >>chng 03 feb 25, this should not be counted here since 00133 * total cooling (less heating) was added into cooling stack */ 00134 for( i=0; i < nCORotate; i++ ) 00135 { 00136 thermal.heatl += C12O16Rotate[i].heat; 00137 thermal.heatl += C13O16Rotate[i].heat; 00138 } 00139 # endif 00140 00141 /* line heating added in following, only here */ 00142 thermal.heating[0][22] = thermal.heatl; 00143 { 00144 enum {DEBUG_LOC=false}; 00145 if( DEBUG_LOC && thermal.heatl/thermal.ctot > 0.1 ) 00146 { 00147 double thresh = 0.1; 00148 fprintf(ioQQQ," all heating lines > %.4f of heatl printed next \n", 00149 thresh ); 00150 for( i=0; i < nWindLine; i++ ) 00151 { 00152 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00153 { 00154 if( TauLine2[i].Coll.heat/thermal.heatl > thresh ) 00155 DumpLine( &TauLine2[i] ); 00156 } 00157 } 00158 00159 for( i=1; i <= nLevel1; i++ ) 00160 { 00161 if( TauLines[i].Coll.heat/thermal.heatl > thresh ) 00162 DumpLine( &TauLines[i] ); 00163 } 00164 00165 /*Atomic & Molecular Lines Chianti & Leiden lines*/ 00166 for( i=0; i <linesAdded2; i++) 00167 { 00168 if( atmolEmis[i].tran->Coll.heat/thermal.heatl > thresh ) 00169 DumpLine( atmolEmis[i].tran ); 00170 } 00171 } 00172 } 00173 00174 /*begin sanity check */ 00175 if( *total <= 0. ) 00176 { 00177 fprintf( ioQQQ, " CoolSum finds cooling <= 0%10.2e\n", 00178 *total ); 00179 } 00180 if( thermal.heatl/thermal.ctot < -1e-15 ) 00181 { 00182 fprintf( ioQQQ, " CoolSum finds negative heating %10.2e %10.2e\n", 00183 thermal.heatl, thermal.ctot ); 00184 } 00185 /*end sanity check */ 00186 00187 thermal.lgCoolEvalOK = false; 00188 return; 00189 }