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 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
00055 hcmap.lgMapDone = true;
00056
00057
00058
00059 PresTotCurrent();
00060
00061
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
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
00110
00111
00112 fac = 1.5;
00113 tmin = torg/fac;
00114 tmax = torg*fac;
00115
00116
00117
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
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
00130 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
00131 double TeNew;
00132 if( thermal.lgTeHigh )
00133 {
00134
00135 factor = 1./factor;
00136
00137 TeNew = (MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)/factor);
00138 }
00139
00140 else
00141 {
00142
00143 TeNew = (tlowst/factor);
00144 }
00145 TempChange(TeNew , false);
00146 }
00147
00148 else
00149 {
00150
00151 fprintf( ioQQQ, " PUNT called with insane argument,=%4.4s\n",
00152 chType );
00153 cdEXIT(EXIT_FAILURE);
00154 }
00155
00156
00157 if( hcmap.nMapAlloc==0 )
00158 {
00159
00160 hcmap.nMapAlloc = hcmap.nMapStep+4;
00161
00162
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
00196 PresTotCurrent();
00197
00198
00199
00200 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00201 {
00202 if( dense.lgElmtOn[nelem] )
00203 {
00204 dense.IonLow[nelem] = 0;
00205
00206
00207
00208
00209
00210
00211 dense.IonHigh[nelem] = nelem + 1;
00212 }
00213 else
00214 {
00215
00216 dense.IonLow[nelem] = -1;
00217 dense.IonHigh[nelem] = -1;
00218 }
00219 }
00220
00221
00222 conv.lgSearch = true;
00223
00224 if( trace.nTrConvg )
00225 fprintf(ioQQQ, " MAP new temp %.4e ===============================================\n",
00226 phycon.te );
00227
00228
00229
00230 conv.nTotalIoniz = 0;
00231
00232
00233 co.iteration_co = -2;
00234
00235
00236 if( ConvEdenIoniz() )
00237 lgAbort = true;
00238
00239
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
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
00278
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
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
00307
00308 fprintf(io, PrintEfmt("%11.4e", phycon.te ) );
00309 fprintf(io, PrintEfmt("%11.4e", thermal.htot ) );
00310 fprintf(io," [%2ld][%2ld]%6.3f",
00311 ksav, jsav,
00312 hfrac);
00313 fprintf(io, PrintEfmt("%11.4e", thermal.ctot ) );
00314 fprintf(io," %s %.1f%c%6.3f",
00315 chLabel ,
00316 wl,
00317 units,
00318 cfrac );
00319 fprintf(io, PrintEfmt(" %10.2e", thermal.dHeatdT ) );
00320 fprintf(io, PrintEfmt("%11.2e", thermal.dCooldT ) );
00321 fprintf(io, PrintEfmt("%11.4e", dense.eden ) );
00322 fprintf(io, 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
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
00356
00357 hcmap.lgMapOK = true;
00358
00359 for( i=2; i< hcmap.nmap-2; ++i )
00360 {
00361 realnum s1,s2,s3;
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
00368
00369
00370
00371
00372
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
00398
00399
00400
00401
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 }