00001
00002
00003
00004
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 t_hcmap hcmap;
00022
00023 void map_do(
00024 FILE *io,
00025 const char *chType)
00026 {
00027 char chLabel[NCOLNT_LAB_LEN+1];
00028 char units;
00029 long int i,
00030 ksav,
00031 j,
00032 jsav,
00033 k,
00034 nelem;
00035 realnum wl;
00036 double cfrac,
00037 ch,
00038 fac,
00039 factor,
00040 hfrac,
00041 oldch,
00042 ratio,
00043 strhet,
00044 strong,
00045 t1,
00046 tlowst,
00047 tmax,
00048 tmin,
00049 torg;
00050
00051 DEBUG_ENTRY( "map_do()" );
00052
00053 t1 = phycon.te;
00054 torg = phycon.te;
00055 hcmap.lgMapOK = true;
00056
00057 hcmap.lgMapDone = true;
00058
00059
00060
00061 PresTotCurrent();
00062
00063
00064 if( called.lgTalk )
00065 {
00066 fprintf( io, " Cloudy punts, Te=%10.3e HTOT=%10.3e CTOT=%10.3e nzone=%4ld\n",
00067 phycon.te, thermal.htot, thermal.ctot, nzone );
00068 fprintf( io, " COOLNG array is\n" );
00069
00070 if( thermal.ctot > 0. )
00071 {
00072 coolpr(io, "ZERO",1,0.,"ZERO");
00073 for( i=0; i < thermal.ncltot; i++ )
00074 {
00075 ratio = thermal.cooling[i]/thermal.ctot;
00076 if( ratio>EPS )
00077 {
00078 coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
00079 ratio,"DOIT");
00080 }
00081 }
00082
00083 coolpr(io, "DONE",1,0.,"DONE");
00084 fprintf( io, " Line heating array follows\n" );
00085 coolpr(io, "ZERO",1,0.,"ZERO");
00086
00087 for( i=0; i < thermal.ncltot; i++ )
00088 {
00089 ratio = thermal.heatnt[i]/thermal.ctot;
00090 if( ratio>EPS )
00091 {
00092 coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
00093 ratio,"DOIT");
00094 }
00095 }
00096
00097 coolpr(io,"DONE",1,0.,"DONE");
00098 }
00099 }
00100
00101
00102 if( called.lgTalk )
00103 {
00104 fprintf( io, " map of heating, cooling, vs temp, follows.\n");
00105 fprintf( io,
00106 " Te Heat--------------------> Cool---------------------> dH/dT dC/DT Ne NH HII Helium \n" );
00107 }
00108
00109 if( strcmp(chType,"punt") == 0 )
00110 {
00111
00112
00113
00114 fac = 1.5;
00115 tmin = torg/fac;
00116 tmax = torg*fac;
00117
00118
00119
00120 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
00121 TempChange(tmin , false);
00122 }
00123
00124 else if( strcmp(chType," map") == 0 )
00125 {
00126
00127 tlowst = MAX2(hcmap.RangeMap[0],phycon.TEMP_LIMIT_LOW);
00128 tmin = tlowst*0.998;
00129 tmax = MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)*1.002;
00130
00131
00132 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
00133 double TeNew;
00134 if( thermal.lgTeHigh )
00135 {
00136
00137 factor = 1./factor;
00138
00139 TeNew = (MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)/factor);
00140 }
00141
00142 else
00143 {
00144
00145 TeNew = (tlowst/factor);
00146 }
00147 TempChange(TeNew , false);
00148 }
00149
00150 else
00151 {
00152
00153 fprintf( ioQQQ, " PUNT called with insane argument,=%4.4s\n",
00154 chType );
00155 cdEXIT(EXIT_FAILURE);
00156 }
00157
00158
00159 if( hcmap.nMapAlloc==0 )
00160 {
00161
00162 hcmap.nMapAlloc = hcmap.nMapStep+4;
00163
00164
00165 hcmap.temap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
00166 if( hcmap.temap == NULL )
00167 {
00168 printf( " not enough memory to allocate hcmap.temap in map_do\n" );
00169 cdEXIT(EXIT_FAILURE);
00170 }
00171 hcmap.cmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
00172 if( hcmap.cmap == NULL )
00173 {
00174 printf( " not enough memory to allocate hcmap.cmap in map_do\n" );
00175 cdEXIT(EXIT_FAILURE);
00176 }
00177 hcmap.hmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
00178 if( hcmap.hmap == NULL )
00179 {
00180 printf( " not enough memory to allocate hcmap.hmap in map_do\n" );
00181 cdEXIT(EXIT_FAILURE);
00182 }
00183 }
00184
00185 thermal.lgCNegChk = false;
00186 hcmap.nmap = 0;
00187 oldch = 0.;
00188 TempChange(phycon.te *factor , true);
00189 if( trace.nTrConvg )
00190 fprintf(ioQQQ, " MAP called temp range %.4e %.4e in %li stops ===============================================\n",
00191 tmin,
00192 tmax,
00193 hcmap.nmap);
00194
00195 while( (double)phycon.te < tmax*0.999 && (double)phycon.te > tmin*1.001 )
00196 {
00197
00198 PresTotCurrent();
00199
00200
00201
00202 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00203 {
00204 if( dense.lgElmtOn[nelem] )
00205 {
00206 dense.IonLow[nelem] = 0;
00207
00208
00209
00210
00211
00212
00213 dense.IonHigh[nelem] = nelem + 1;
00214 }
00215 else
00216 {
00217
00218 dense.IonLow[nelem] = -1;
00219 dense.IonHigh[nelem] = -1;
00220 }
00221 }
00222
00223
00224 conv.lgSearch = true;
00225
00226 if( trace.nTrConvg )
00227 fprintf(ioQQQ, " MAP new temp %.4e ===============================================\n",
00228 phycon.te );
00229
00230
00231
00232 conv.nTotalIoniz = 0;
00233
00234
00235 if( ConvEdenIoniz() )
00236 lgAbort = true;
00237
00238
00239 hcmap.temap[hcmap.nmap] = (realnum)phycon.te;
00240 hcmap.cmap[hcmap.nmap] = (realnum)thermal.ctot;
00241 hcmap.hmap[hcmap.nmap] = (realnum)thermal.htot;
00242
00243 wl = 0.f;
00244 strong = 0.;
00245
00246 for( j=0; j < thermal.ncltot; j++ )
00247 {
00248 if( thermal.cooling[j] > strong )
00249 {
00250 strcpy( chLabel, thermal.chClntLab[j] );
00251 strong = thermal.cooling[j];
00252 wl = thermal.collam[j];
00253 }
00254 }
00255
00256 cfrac = strong/thermal.ctot;
00257 strhet = 0.;
00258
00259 ksav = -INT_MAX;
00260 jsav = -INT_MAX;
00261
00262 for( k=0; k < LIMELM; k++ )
00263 {
00264 for( j=0; j < LIMELM; j++ )
00265 {
00266 if( thermal.heating[k][j] > strhet )
00267 {
00268 strhet = thermal.heating[k][j];
00269 jsav = j;
00270 ksav = k;
00271 }
00272 }
00273 }
00274
00275 ch = thermal.ctot - thermal.htot;
00276
00277
00278 if( oldch/ch < 0. && called.lgTalk )
00279 {
00280 fprintf( io, " ----------------------------------------------- Probable thermal solution here. --------------------------------------------\n" );
00281 }
00282
00283 oldch = ch;
00284 hfrac = strhet/thermal.htot;
00285 if( called.lgTalk )
00286 {
00287
00288 if( wl < 100000.f )
00289 {
00290 units = 'A';
00291 }
00292
00293 else
00294 {
00295 wl /= 10000.f;
00296 units = 'm';
00297 }
00298
00299 if( trace.lgTrace )
00300 {
00301 fprintf( io, "TRACE: te, htot, ctot%11.3e%11.3e%11.3e\n",
00302 phycon.te, thermal.htot, thermal.ctot );
00303 }
00304
00305
00306
00307 fprintf(io, PrintEfmt("%11.4e", phycon.te ) );
00308 fprintf(io, PrintEfmt("%11.4e", thermal.htot ) );
00309 fprintf(io," [%2ld][%2ld]%6.3f",
00310 ksav, jsav,
00311 hfrac);
00312 fprintf(io, PrintEfmt("%11.4e", thermal.ctot ) );
00313 fprintf(io," %s %.1f%c%6.3f",
00314 chLabel ,
00315 wl,
00316 units,
00317 cfrac );
00318 fprintf(io, PrintEfmt(" %10.2e", thermal.dHeatdT ) );
00319 fprintf(io, PrintEfmt("%11.2e", thermal.dCooldT ) );
00320 fprintf(io, PrintEfmt("%11.4e", dense.eden ) );
00321 fprintf(io, PrintEfmt("%11.4e", dense.gas_phase[ipHYDROGEN] ) );
00322 if( dense.lgElmtOn[ipHELIUM] )
00323 {
00324 fprintf(io,"%6.2f%6.2f%6.2f%6.2f\n",
00325 log10(MAX2(1e-9,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])),
00326 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][0]/dense.gas_phase[ipHELIUM])),
00327 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM])),
00328 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][2]/dense.gas_phase[ipHELIUM])) );
00329 }
00330 fflush(io);
00331 }
00332
00333 TempChange(phycon.te*factor , true);
00334
00335 hcmap.nmap = MIN2(hcmap.nMapAlloc,hcmap.nmap+1);
00336
00337 {
00338 enum {DEBUG_LOC=false};
00339 if( DEBUG_LOC )
00340 {
00341 static int kount = 0;
00342 factor = 1.;
00343 TempChange(8674900. , true);
00344 ++kount;
00345 if( kount >=100 )
00346 {
00347 fprintf(ioQQQ," exiting in map_do\n");
00348 break;
00349 }
00350 }
00351 }
00352 }
00353
00354
00355
00356 hcmap.lgMapOK = true;
00357
00358 for( i=2; i< hcmap.nmap-2; ++i )
00359 {
00360 realnum s1,s2,s3;
00361 s1 = hcmap.cmap[i-2] - hcmap.cmap[i-1];
00362 s2 = hcmap.cmap[i-1] - hcmap.cmap[i];
00363 s3 = hcmap.cmap[i] - hcmap.cmap[i+1];
00364 if( s1*s3 > 0. && s2*s3 < 0. )
00365 {
00366
00367
00368
00369
00370
00371
00372 fprintf( io,
00373 " cooling curve had double inflection at T=%.2e. ",
00374 hcmap.temap[i]);
00375 fprintf( io, " Slopes were %.2e %.2e %.2e", s1, s2, s3);
00376 if( fabs(s2)/hcmap.cmap[i] > 0.05 )
00377 {
00378 fprintf( io,
00379 " error large, (rel slope of %.2e).\n",
00380 s2 / hcmap.cmap[i]);
00381 hcmap.lgMapOK = false;
00382 }
00383 else
00384 {
00385 fprintf( io,
00386 " error is small, (rel slope of %.2e).\n",
00387 s2 / hcmap.cmap[i]);
00388 }
00389 }
00390
00391 s1 = hcmap.hmap[i-2] - hcmap.hmap[i-1];
00392 s2 = hcmap.hmap[i-1] - hcmap.hmap[i];
00393 s3 = hcmap.hmap[i] - hcmap.hmap[i+1];
00394 if( s1*s3 > 0. && s2*s3 < 0. )
00395 {
00396
00397
00398
00399
00400
00401 fprintf( io,
00402 " heating curve had double inflection at T=%.2e.\n",
00403 hcmap.temap[i] );
00404 hcmap.lgMapOK = false;
00405 }
00406 }
00407
00408 thermal.lgCNegChk = true;
00409 TempChange(t1 , false);
00410 return;
00411 }