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