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 /*HeatPunch punch contributors to local heating, with punch heat command, called by punch_do */ 00004 #include "cddefines.h" 00005 #include "thermal.h" 00006 #include "radius.h" 00007 #include "conv.h" 00008 #include "lines_service.h" 00009 #include "dense.h" 00010 #include "taulines.h" 00011 #include "phycon.h" 00012 #include "elementnames.h" 00013 #include "dynamics.h" 00014 #include "punch.h" 00015 00016 /* limit for number of heat agents that are saved */ 00017 /* limit to number to print */ 00018 static const int IPRINT= 9; 00019 00020 /*HeatPunch punch contributors to local heating, with punch heat command, 00021 * called by punch_do */ 00022 void HeatPunch(FILE* io) 00023 { 00024 char **chLabel, 00025 chLbl[11]; 00026 bool lgHeatLine; 00027 int nFail; 00028 long int i, 00029 ipStrong, 00030 ipnt, 00031 *ipOrdered, 00032 *ipsave, 00033 j, 00034 *jpsave, 00035 k, 00036 level; 00037 double CS, 00038 ColHeat, 00039 EscP, 00040 Pump, 00041 Strong, 00042 TauIn, 00043 cool_total, 00044 heat_total; 00045 realnum *save; 00046 00047 DEBUG_ENTRY( "HeatPunch()" ); 00048 00049 save = (realnum *) CALLOC(LIMELM*LIMELM, sizeof(realnum)); 00050 ipsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int)); 00051 jpsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int)); 00052 ipOrdered = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int)); 00053 chLabel = (char **) CALLOC(LIMELM*LIMELM, sizeof(char *)); 00054 00055 for( i=0; i<LIMELM*LIMELM; ++i ) 00056 { 00057 ipsave[i] = INT_MIN; 00058 jpsave[i] = INT_MIN; 00059 save[i] = -FLT_MAX; 00060 chLabel[i] = (char *) CALLOC(10, sizeof(char)); 00061 } 00062 00063 cool_total = thermal.ctot; 00064 heat_total = thermal.htot; 00065 00066 /* >>chng 06 mar 17, comment out following block and replace with this 00067 * removing dynamics heating & cooling and report only physical 00068 * heating and cooling 00069 * NB the heating and cooling as punched no longer need be 00070 * equal for a converged model */ 00071 cool_total -= dynamics.Cool; 00072 heat_total -= dynamics.Heat; 00073 # if 0 00074 if(dynamics.Cool > dynamics.Heat) 00075 { 00076 cool_total -= dynamics.Heat; 00077 heat_total -= dynamics.Heat; 00078 } 00079 else 00080 { 00081 cool_total -= dynamics.Cool; 00082 heat_total -= dynamics.Cool; 00083 } 00084 # endif 00085 00086 ipnt = 0; 00087 00088 /* heat sources are saved in a 2d square array */ 00089 /* WeakHeatCool set with 'set weakheatcool' command 00090 * default is 0.05 */ 00091 for( i=0; i < LIMELM; i++ ) 00092 { 00093 for( j=0; j < LIMELM; j++ ) 00094 { 00095 if( thermal.heating[i][j]/heat_total > SMALLFLOAT ) 00096 { 00097 ipsave[ipnt] = i; 00098 jpsave[ipnt] = j; 00099 save[ipnt] = (realnum)(thermal.heating[i][j]/heat_total); 00100 ipnt++; 00101 } 00102 } 00103 } 00104 00105 /* now check for possible line heating not counted in 1,23 00106 * there should not be any significant heating source here 00107 * since they would not be counted in derivative correctly */ 00108 for( i=0; i < thermal.ncltot; i++ ) 00109 { 00110 if( thermal.heatnt[i]/heat_total > punch.WeakHeatCool ) 00111 { 00112 realnum awl; 00113 awl = thermal.collam[i]; 00114 /* force to save wavelength convention as printout */ 00115 if( awl > 100000 ) 00116 awl /= 10000; 00117 fprintf( io, " Negative coolant was %s %.2f %.2e\n", 00118 thermal.chClntLab[i], awl, thermal.heatnt[i]/heat_total ); 00119 } 00120 } 00121 00122 if( !conv.lgConvTemp ) 00123 { 00124 fprintf( io, "#>>>> Temperature not converged.\n" ); 00125 } 00126 else if( !conv.lgConvEden ) 00127 { 00128 fprintf( io, "#>>>> Electron density not converged.\n" ); 00129 } 00130 else if( !conv.lgConvIoniz ) 00131 { 00132 fprintf( io, "#>>>> Ionization not converged.\n" ); 00133 } 00134 else if( !conv.lgConvPres ) 00135 { 00136 fprintf( io, "#>>>> Pressure not converged.\n" ); 00137 } 00138 00139 /* this is mainly to keep the compiler from flagging possible paths 00140 * with j not being set */ 00141 i = INT_MIN; 00142 j = INT_MIN; 00143 /* following loop tries to identify strongest agents and turn to labels */ 00144 for( k=0; k < ipnt; k++ ) 00145 { 00146 /* generate labels that make sense in printout 00147 * if not identified with a specific agent, will print indices as [i][j], 00148 * heating is thermal.heating[i][j] */ 00149 i = ipsave[k]; 00150 j = jpsave[k]; 00151 /* i >= j indicates agent is one of the first LIMELM elements */ 00152 if( i >= j ) 00153 { 00154 if( dense.xIonDense[i][j] == 0. && thermal.heating[i][j]>SMALLFLOAT ) 00155 fprintf(ioQQQ,"DISASTER assert about to be thrown - search for hit it\n"); 00156 /*fprintf(ioQQQ,"DEBUG %li %li %.2e %.2e\n", i , j , 00157 dense.xIonDense[i][j], 00158 thermal.heating[i][j]);*/ 00159 ASSERT( dense.xIonDense[i][j] > 0. || thermal.heating[i][j]<SMALLFLOAT ); 00160 /* this is case of photoionization of atom or ion */ 00161 strcpy( chLabel[k], elementnames.chElementSym[i] ); 00162 strcat( chLabel[k], elementnames.chIonStage[j] ); 00163 } 00164 /* notice that in test i and j are swapped from order in heating array */ 00165 else if( i == 0 && j == 1 ) 00166 { 00167 /* photoionization from all excited states of Hydrogenic species, 00168 * heating[0][1] */ 00169 strcpy( chLabel[k], "Hn=2" ); 00170 } 00171 else if( i == 0 && j == 3 ) 00172 { 00173 /* collisional ionization of all iso-seq from all levels, 00174 * heating[0][3] */ 00175 strcpy( chLabel[k], "Hion" ); 00176 } 00177 else if( i == 0 && j == 7 ) 00178 { 00179 /* UTA ionization heating[0][7] */ 00180 strcpy( chLabel[k], " UTA" ); 00181 } 00182 else if( i == 0 && j == 8 ) 00183 { 00184 /* thermal.heating[0][8] is heating due to collisions within 00185 * X of H2, code var hmi.HeatH2Dexc_used, hmi.HeatH2Dexc_BigH2, 00186 * hmi.HeatH2Dexc_TH85 */ 00187 strcpy( chLabel[k], "H2vH" ); 00188 } 00189 else if( i == 0 && j == 17 ) 00190 { 00191 /* thermal.heating[0][17] is heating due to photodissociation 00192 * heating of X within H2, 00193 * code var hmi.HeatH2Dish_used */ 00194 strcpy( chLabel[k], "H2dH" ); 00195 } 00196 else if( i == 0 && j == 9 ) 00197 { 00198 /* CO dissociation, co.CODissHeat, heating[0][9] */ 00199 strcpy( chLabel[k], "COds" ); 00200 } 00201 else if( i == 0 && j == 20 ) 00202 { 00203 /* extra heat thermal.heating[0][20]*/ 00204 strcpy( chLabel[k], "extH" ); 00205 } 00206 else if( i == 0 && j == 11 ) 00207 { 00208 /* free free heating, heating[0][11] */ 00209 strcpy( chLabel[k], "H FF" ); 00210 } 00211 else if( i == 0 && j == 12 ) 00212 { 00213 /* heating line, that was supposed to cool, heating[0][12] */ 00214 strcpy( chLabel[k], "Clin" ); 00215 } 00216 else if( i == 0 && j == 13 ) 00217 { 00218 /* grain photoionization, heating[0][13] */ 00219 strcpy( chLabel[k], "GrnP" ); 00220 } 00221 else if( i == 0 && j == 14 ) 00222 { 00223 /* grain collisions, heating[0][14] */ 00224 strcpy( chLabel[k], "GrnC" ); 00225 } 00226 else if( i == 0 && j == 15 ) 00227 { 00228 /* H- heating, heating[0][15] */ 00229 strcpy( chLabel[k], "H- " ); 00230 } 00231 else if( i == 0 && j == 16 ) 00232 { 00233 /* H2+ heating, heating[0][16] */ 00234 strcpy( chLabel[k], "H2+ " ); 00235 } 00236 else if( i == 0 && j == 18 ) 00237 { 00238 /* H2 photoionization heating, heating[0][18] */ 00239 strcpy( chLabel[k], "H2ph" ); 00240 } 00241 else if( i == 0 && j == 19 ) 00242 { 00243 /* Compton heating, heating[0][19] */ 00244 strcpy( chLabel[k], "Comp" ); 00245 } 00246 else if( i == 0 && j == 22 ) 00247 { 00248 /* line heating, heating[0][22] */ 00249 strcpy( chLabel[k], "line" ); 00250 } 00251 else if( i == 0 && j == 23 ) 00252 { 00253 /* iso-sequence line heating - all elements together, 00254 * heating[0][23] */ 00255 strcpy( chLabel[k], "Hlin" ); 00256 } 00257 else if( i == 0 && j == 24 ) 00258 { 00259 /* charge transfer heating, heating[0][24] */ 00260 strcpy( chLabel[k], "ChaT" ); 00261 } 00262 else if( i == 1 && j == 3 ) 00263 { 00264 /* helium triplet line heating, heating[1][3] */ 00265 strcpy( chLabel[k], "He3l" ); 00266 } 00267 else if( i == 1 && j == 5 ) 00268 { 00269 /* advective heating, heating[1][5] */ 00270 strcpy( chLabel[k], "adve" ); 00271 } 00272 else if( i == 1 && j == 6 ) 00273 { 00274 /* cosmic ray heating thermal.heating[1][6]*/ 00275 strcpy( chLabel[k], "CR H" ); 00276 } 00277 else if( i == 25 && j == 27 ) 00278 { 00279 /* Fe 2 line heating, heating[25][27] */ 00280 strcpy( chLabel[k], "Fe 2" ); 00281 } 00282 else 00283 { 00284 sprintf( chLabel[k], "[%ld][%ld]" , i , j ); 00285 } 00286 } 00287 00288 /* now sort by decreasing importance */ 00289 /*spsort netlib routine to sort array returning sorted indices */ 00290 spsort( 00291 /* input array to be sorted */ 00292 save, 00293 /* number of values in x */ 00294 ipnt, 00295 /* permutation output array */ 00296 ipOrdered, 00297 /* flag saying what to do - 1 sorts into increasing order, not changing 00298 * the original routine */ 00299 -1, 00300 /* error condition, should be 0 */ 00301 &nFail); 00302 00303 /*>>chng 06 jun 06, change start of punch to give same info as cooling 00304 * as per comment by Yumihiko Tsuzuki */ 00305 /* begin the print out with zone number, total heating and cooling */ 00306 fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e", 00307 radius.depth_mid_zone, 00308 phycon.te, 00309 heat_total, 00310 cool_total ); 00311 00312 /* we only want to print the IPRINT strongest of the group */ 00313 ipnt = MIN2( ipnt , IPRINT ); 00314 00315 for( k=0; k < ipnt; k++ ) 00316 { 00317 int ip = ipOrdered[k]; 00318 i = ipsave[ip]; 00319 j = jpsave[ip]; 00320 ASSERT( i<LIMELM && j<LIMELM ); 00321 if(k > 4 && thermal.heating[i][j]/heat_total < punch.WeakHeatCool) 00322 break; 00323 fprintf( io, "\t%s\t%.7f ", 00324 chLabel[ip], save[ip] ); 00325 } 00326 fprintf( io, " \n" ); 00327 00328 /* a negative pointer in the heating array is probably a problem, 00329 * indicating that some line has become a heat source */ 00330 lgHeatLine = false; 00331 00332 /* check if any lines were major heat sources */ 00333 for( i=0; i < ipnt; i++ ) 00334 { 00335 /* heating[22][0] is line heating - identify line if important */ 00336 if( ipsave[i] == 0 && jpsave[i] == 22 ) 00337 lgHeatLine = true; 00338 } 00339 00340 if( lgHeatLine ) 00341 { 00342 /* a line was a major heat source - identify it */ 00343 FndLineHt(&level,&ipStrong,&Strong); 00344 if( Strong/heat_total > 0.005 ) 00345 { 00346 if( level == 1 ) 00347 { 00348 strcpy( chLbl, chLineLbl(&TauLines[ipStrong] ) ); 00349 TauIn = TauLines[ipStrong].Emis->TauIn; 00350 Pump = TauLines[ipStrong].Emis->pump; 00351 EscP = TauLines[ipStrong].Emis->Pesc; 00352 CS = TauLines[ipStrong].Coll.cs; 00353 /* ratio of line to total heating */ 00354 ColHeat = TauLines[ipStrong].Coll.heat/ 00355 heat_total; 00356 } 00357 else if( level == 2 ) 00358 { 00359 strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]) ); 00360 TauIn = TauLine2[ipStrong].Emis->TauIn; 00361 Pump = TauLine2[ipStrong].Emis->pump; 00362 EscP = TauLine2[ipStrong].Emis->Pesc; 00363 CS = TauLine2[ipStrong].Coll.cs; 00364 ColHeat = TauLine2[ipStrong].Coll.heat/ 00365 heat_total; 00366 } 00367 else if( level == 3 ) 00368 /* C12O16 */ 00369 { 00370 strcpy( chLbl, chLineLbl(&C12O16Rotate[ipStrong]) ); 00371 TauIn = C12O16Rotate[ipStrong].Emis->TauIn; 00372 Pump = C12O16Rotate[ipStrong].Emis->pump; 00373 EscP = C12O16Rotate[ipStrong].Emis->Pesc; 00374 CS = C12O16Rotate[ipStrong].Coll.cs; 00375 ColHeat = C12O16Rotate[ipStrong].Coll.heat/ 00376 heat_total; 00377 } 00378 else if( level == 4 ) 00379 /* C13O16 */ 00380 { 00381 strcpy( chLbl, chLineLbl(&C13O16Rotate[ipStrong]) ); 00382 TauIn = C13O16Rotate[ipStrong].Emis->TauIn; 00383 Pump = C13O16Rotate[ipStrong].Emis->pump; 00384 EscP = C13O16Rotate[ipStrong].Emis->Pesc; 00385 CS = C13O16Rotate[ipStrong].Coll.cs; 00386 ColHeat = C13O16Rotate[ipStrong].Coll.heat/ 00387 heat_total; 00388 } 00389 else if( level == 5 ) 00390 /* hyperfine levels */ 00391 { 00392 strcpy( chLbl, chLineLbl(&HFLines[ipStrong]) ); 00393 TauIn = HFLines[ipStrong].Emis->TauIn; 00394 Pump = HFLines[ipStrong].Emis->pump; 00395 EscP = HFLines[ipStrong].Emis->Pesc; 00396 CS = HFLines[ipStrong].Coll.cs; 00397 ColHeat = HFLines[ipStrong].Coll.heat/ 00398 heat_total; 00399 } 00400 else 00401 { 00402 fprintf( ioQQQ, " HeatPunch insane level%4ld\n", 00403 level ); 00404 cdEXIT(EXIT_FAILURE); 00405 } 00406 fprintf( io, " LHeat lv%2ld %10.10s TIn%10.2e Pmp%9.1e EscP%9.1e CS%9.1e Hlin/tot%10.2e\n", 00407 level, chLbl, TauIn, Pump, EscP, CS, ColHeat ); 00408 } 00409 } 00410 for( i=0; i<LIMELM*LIMELM; ++i ) 00411 { 00412 free(chLabel[i]); 00413 } 00414 00415 free(chLabel); 00416 free(ipOrdered); 00417 free(jpsave); 00418 free(ipsave); 00419 free(save); 00420 return; 00421 }