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 /*map_do produce map of heating-cooling space for specified zone, called as result of 00004 * map command */ 00005 #include "cddefines.h" 00006 #include "thermal.h" 00007 #include "cooling.h" 00008 #include "called.h" 00009 #include "dense.h" 00010 #include "mole.h" 00011 #include "phycon.h" 00012 #include "trace.h" 00013 #include "pressure.h" 00014 #include "conv.h" 00015 #include "hcmap.h" 00016 #ifdef EPS 00017 # undef EPS 00018 #endif 00019 #define EPS 0.005 00020 00021 void map_do( 00022 FILE *io, 00023 const char *chType) 00024 { 00025 char chLabel[NCOLNT_LAB_LEN+1]; 00026 char units; 00027 long int i, 00028 ksav, 00029 j, 00030 jsav, 00031 k, 00032 nelem; 00033 realnum wl; 00034 double cfrac, 00035 ch, 00036 fac, 00037 factor, 00038 hfrac, 00039 oldch, 00040 ratio, 00041 strhet, 00042 strong, 00043 t1, 00044 tlowst, 00045 tmax, 00046 tmin, 00047 torg; 00048 00049 DEBUG_ENTRY( "map_do()" ); 00050 00051 t1 = phycon.te; 00052 torg = phycon.te; 00053 hcmap.lgMapOK = true; 00054 /* flag indicating that we have computed a map */ 00055 hcmap.lgMapDone = true; 00056 00057 /* make sure pressure has been evaluated */ 00058 /* this sets values of pressure.PresTotlCurr */ 00059 PresTotCurrent(); 00060 00061 /* print out all coolants if all else fails */ 00062 if( called.lgTalk ) 00063 { 00064 fprintf( io, " Cloudy punts, Te=%10.3e HTOT=%10.3e CTOT=%10.3e nzone=%4ld\n", 00065 phycon.te, thermal.htot, thermal.ctot, nzone ); 00066 fprintf( io, " COOLNG array is\n" ); 00067 00068 if( thermal.ctot > 0. ) 00069 { 00070 coolpr(io, "ZERO",1,0.,"ZERO"); 00071 for( i=0; i < thermal.ncltot; i++ ) 00072 { 00073 ratio = thermal.cooling[i]/thermal.ctot; 00074 if( ratio>EPS ) 00075 { 00076 coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i], 00077 ratio,"DOIT"); 00078 } 00079 } 00080 00081 coolpr(io, "DONE",1,0.,"DONE"); 00082 fprintf( io, " Line heating array follows\n" ); 00083 coolpr(io, "ZERO",1,0.,"ZERO"); 00084 00085 for( i=0; i < thermal.ncltot; i++ ) 00086 { 00087 ratio = thermal.heatnt[i]/thermal.ctot; 00088 if( ratio>EPS ) 00089 { 00090 coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i], 00091 ratio,"DOIT"); 00092 } 00093 } 00094 00095 coolpr(io,"DONE",1,0.,"DONE"); 00096 } 00097 } 00098 00099 /* map out te-ionization-cooling space before punching out. */ 00100 if( called.lgTalk ) 00101 { 00102 fprintf( io, " map of heating, cooling, vs temp, follows.\n"); 00103 fprintf( io, 00104 " Te Heat--------------------> Cool---------------------> dH/dT dC/DT Ne NH HII Helium \n" ); 00105 } 00106 00107 if( strcmp(chType,"punt") == 0 ) 00108 { 00109 /* this is the original use of punt, we are punting 00110 * only map within factor of two of final temperature 00111 * fac will the range to either side of punted temperature */ 00112 fac = 1.5; 00113 tmin = torg/fac; 00114 tmax = torg*fac; 00115 00116 /* we want about 20 steps between high and low temperature 00117 * default of nMapStep is 20, set with set nmaps command */ 00118 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep)); 00119 TempChange(tmin , false); 00120 } 00121 00122 else if( strcmp(chType," map") == 0 ) 00123 { 00124 /* create some sort of map of heating-cooling */ 00125 tlowst = MAX2(hcmap.RangeMap[0],phycon.TEMP_LIMIT_LOW); 00126 tmin = tlowst*0.998; 00127 tmax = MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)*1.002; 00128 00129 /* we want about nMapStep (=20) steps between high and low temperature */ 00130 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep)); 00131 double TeNew; 00132 if( thermal.lgTeHigh ) 00133 { 00134 /* high te */ 00135 factor = 1./factor; 00136 /* TeHighest is highest possible temperature, 1E10 */ 00137 TeNew = (MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)/factor); 00138 } 00139 00140 else 00141 { 00142 /* low te */ 00143 TeNew = (tlowst/factor); 00144 } 00145 TempChange(TeNew , false); 00146 } 00147 00148 else 00149 { 00150 /* don't know what to do */ 00151 fprintf( ioQQQ, " PUNT called with insane argument,=%4.4s\n", 00152 chType ); 00153 cdEXIT(EXIT_FAILURE); 00154 } 00155 00156 /* now allocate space for te, c, h vectors in map, if not already done */ 00157 if( hcmap.nMapAlloc==0 ) 00158 { 00159 /* space not allocated, do so now */ 00160 hcmap.nMapAlloc = hcmap.nMapStep+4; 00161 00162 /* now make the space */ 00163 hcmap.temap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) ); 00164 if( hcmap.temap == NULL ) 00165 { 00166 printf( " not enough memory to allocate hcmap.temap in map_do\n" ); 00167 cdEXIT(EXIT_FAILURE); 00168 } 00169 hcmap.cmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) ); 00170 if( hcmap.cmap == NULL ) 00171 { 00172 printf( " not enough memory to allocate hcmap.cmap in map_do\n" ); 00173 cdEXIT(EXIT_FAILURE); 00174 } 00175 hcmap.hmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) ); 00176 if( hcmap.hmap == NULL ) 00177 { 00178 printf( " not enough memory to allocate hcmap.hmap in map_do\n" ); 00179 cdEXIT(EXIT_FAILURE); 00180 } 00181 } 00182 00183 thermal.lgCNegChk = false; 00184 hcmap.nmap = 0; 00185 oldch = 0.; 00186 TempChange(phycon.te *factor , true); 00187 if( trace.nTrConvg ) 00188 fprintf(ioQQQ, " MAP called temp range %.4e %.4e in %li stops ===============================================\n", 00189 tmin, 00190 tmax, 00191 hcmap.nmap); 00192 00193 while( (double)phycon.te < tmax*0.999 && (double)phycon.te > tmin*1.001 ) 00194 { 00195 /* this sets values of pressure.PresTotlCurr */ 00196 PresTotCurrent(); 00197 00198 /* must reset upper and lower bounds for ionization distributions */ 00199 /* fix number of stages of ionization */ 00200 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00201 { 00202 if( dense.lgElmtOn[nelem] ) 00203 { 00204 dense.IonLow[nelem] = 0; 00205 /* 00206 * IonHigh[n] is the highest stage of ionization present 00207 * the IonHigh array index is on the C scake, so [0] is hydrogen 00208 * the value is also on the C scale, so element [nelem] can range 00209 * from 0 to nelem+1 00210 */ 00211 dense.IonHigh[nelem] = nelem + 1; 00212 } 00213 else 00214 { 00215 /* this element is turned off, make stages impossible */ 00216 dense.IonLow[nelem] = -1; 00217 dense.IonHigh[nelem] = -1; 00218 } 00219 } 00220 00221 /* this turns on constant reevaluation of everything */ 00222 conv.lgSearch = true; 00223 00224 if( trace.nTrConvg ) 00225 fprintf(ioQQQ, " MAP new temp %.4e ===============================================\n", 00226 phycon.te ); 00227 00228 /* this counts how many times ionize is called in this model after startr, 00229 * and is flag used by ionize to understand it is being called the first time*/ 00230 conv.nTotalIoniz = 0; 00231 00232 /* this will force reinitialization of co network */ 00233 co.iteration_co = -2; 00234 00235 /* now get ionization solution for this temperature */ 00236 if( ConvEdenIoniz() ) 00237 lgAbort = true; 00238 00239 /* save results for later prints */ 00240 hcmap.temap[hcmap.nmap] = (realnum)phycon.te; 00241 hcmap.cmap[hcmap.nmap] = (realnum)thermal.ctot; 00242 hcmap.hmap[hcmap.nmap] = (realnum)thermal.htot; 00243 00244 wl = 0.f; 00245 strong = 0.; 00246 00247 for( j=0; j < thermal.ncltot; j++ ) 00248 { 00249 if( thermal.cooling[j] > strong ) 00250 { 00251 strcpy( chLabel, thermal.chClntLab[j] ); 00252 strong = thermal.cooling[j]; 00253 wl = thermal.collam[j]; 00254 } 00255 } 00256 00257 cfrac = strong/thermal.ctot; 00258 strhet = 0.; 00259 /* these will be reset in following loop*/ 00260 ksav = -INT_MAX; 00261 jsav = -INT_MAX; 00262 00263 for( k=0; k < LIMELM; k++ ) 00264 { 00265 for( j=0; j < LIMELM; j++ ) 00266 { 00267 if( thermal.heating[k][j] > strhet ) 00268 { 00269 strhet = thermal.heating[k][j]; 00270 jsav = j; 00271 ksav = k; 00272 } 00273 } 00274 } 00275 00276 ch = thermal.ctot - thermal.htot; 00277 /* use ratio to check for change of sign since product 00278 * can underflow at low densities */ 00279 if( oldch/ch < 0. && called.lgTalk ) 00280 { 00281 fprintf( io, " ----------------------------------------------- Probable thermal solution here. --------------------------------------------\n" ); 00282 } 00283 00284 oldch = ch; 00285 hfrac = strhet/thermal.htot; 00286 if( called.lgTalk ) 00287 { 00288 /* convert to micros if IR transition */ 00289 if( wl < 100000.f ) 00290 { 00291 units = 'A'; 00292 } 00293 00294 else 00295 { 00296 wl /= 10000.f; 00297 units = 'm'; 00298 } 00299 00300 if( trace.lgTrace ) 00301 { 00302 fprintf( io, "TRACE: te, htot, ctot%11.3e%11.3e%11.3e\n", 00303 phycon.te, thermal.htot, thermal.ctot ); 00304 } 00305 00306 /*fprintf( io, 00307 "%10.4e%11.4e%4ld%4ld%6.3f%11.4e %4.4s %4ld%c%6.3f%10.2e%11.4e%11.4e%6.2f%6.2f%6.2f%6.2f\n",;*/ 00308 fprintf(io,"%s", PrintEfmt("%11.4e", phycon.te ) ); 00309 fprintf(io,"%s", PrintEfmt("%11.4e", thermal.htot ) ); 00310 fprintf(io," [%2ld][%2ld]%6.3f", 00311 ksav, jsav, 00312 hfrac); 00313 fprintf(io,"%s", PrintEfmt("%11.4e", thermal.ctot ) ); 00314 fprintf(io," %s %.1f%c%6.3f", 00315 chLabel , 00316 wl, 00317 units, 00318 cfrac ); 00319 fprintf(io,"%s", PrintEfmt(" %10.2e", thermal.dHeatdT ) ); 00320 fprintf(io,"%s", PrintEfmt("%11.2e", thermal.dCooldT ) ); 00321 fprintf(io,"%s", PrintEfmt("%11.4e", dense.eden ) ); 00322 fprintf(io,"%s", PrintEfmt("%11.4e", dense.gas_phase[ipHYDROGEN] ) ); 00323 if( dense.lgElmtOn[ipHELIUM] ) 00324 { 00325 fprintf(io,"%6.2f%6.2f%6.2f%6.2f\n", 00326 log10(MAX2(1e-9,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])), 00327 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][0]/dense.gas_phase[ipHELIUM])), 00328 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM])), 00329 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][2]/dense.gas_phase[ipHELIUM])) ); 00330 } 00331 fflush(io); 00332 } 00333 00334 TempChange(phycon.te*factor , true); 00335 /* increment nmap but do not exceed nMapAlloc */ 00336 hcmap.nmap = MIN2(hcmap.nMapAlloc,hcmap.nmap+1); 00337 00338 { 00339 enum {DEBUG_LOC=false}; 00340 if( DEBUG_LOC ) 00341 { 00342 static int kount = 0; 00343 factor = 1.; 00344 TempChange(8674900. , true); 00345 ++kount; 00346 if( kount >=100 ) 00347 { 00348 fprintf(ioQQQ," exiting in map_do\n"); 00349 break; 00350 } 00351 } 00352 } 00353 } 00354 00355 /* now check whether sharp inflections occurred, and also find the biggest jump 00356 * in the heating and cooling */ 00357 hcmap.lgMapOK = true; 00358 /* >>chng 02 mar 04, lower bound had been 1, so [i-2] below was negative */ 00359 for( i=2; i< hcmap.nmap-2; ++i ) 00360 { 00361 realnum s1,s2,s3;/* the three slopes we will use */ 00362 s1 = hcmap.cmap[i-2] - hcmap.cmap[i-1]; 00363 s2 = hcmap.cmap[i-1] - hcmap.cmap[i]; 00364 s3 = hcmap.cmap[i] - hcmap.cmap[i+1]; 00365 if( s1*s3 > 0. && s2*s3 < 0. ) 00366 { 00367 /* of the three points, the outer had the same slope 00368 * (their product was positive) but there was an inflection 00369 * between them (the negative product). The data chain looked like 00370 * 2 4 00371 * 1 3 or vice versa, either case is wrong, 00372 * with the logic in the above test, the problem point will aways be s2 */ 00373 fprintf( io, 00374 " cooling curve had double inflection at T=%.2e. ", 00375 hcmap.temap[i]); 00376 fprintf( io, " Slopes were %.2e %.2e %.2e", s1, s2, s3); 00377 if( fabs(s2)/hcmap.cmap[i] > 0.05 ) 00378 { 00379 fprintf( io, 00380 " error large, (rel slope of %.2e).\n", 00381 s2 / hcmap.cmap[i]); 00382 hcmap.lgMapOK = false; 00383 } 00384 else 00385 { 00386 fprintf( io, 00387 " error is small, (rel slope of %.2e).\n", 00388 s2 / hcmap.cmap[i]); 00389 } 00390 } 00391 00392 s1 = hcmap.hmap[i-2] - hcmap.hmap[i-1]; 00393 s2 = hcmap.hmap[i-1] - hcmap.hmap[i]; 00394 s3 = hcmap.hmap[i] - hcmap.hmap[i+1]; 00395 if( s1*s3 > 0. && s2*s3 < 0. ) 00396 { 00397 /* of the three points, the outer had the same slope 00398 * (their product was positive) but there was an inflection 00399 * between them (the negative product). The data chain looked like 00400 * 2 4 00401 * 1 3 or vice versa, either case is wrong */ 00402 fprintf( io, 00403 " heating curve had double inflection at T=%.2e.\n", 00404 hcmap.temap[i] ); 00405 hcmap.lgMapOK = false; 00406 } 00407 } 00408 00409 thermal.lgCNegChk = true; 00410 TempChange(t1 , false); 00411 return; 00412 }