00001
00002
00003
00004
00005
00006
00007
00008 #include "cddefines.h"
00009 #define Z 1.0001
00010 #include "wind.h"
00011 #include "stopcalc.h"
00012 #include "thermal.h"
00013 #include "dynamics.h"
00014 #include "trace.h"
00015 #include "punch.h"
00016 #include "pressure.h"
00017 #include "iso.h"
00018 #include "h2.h"
00019 #include "rfield.h"
00020 #include "dense.h"
00021 #include "hmi.h"
00022 #include "geometry.h"
00023 #include "opacity.h"
00024 #include "ipoint.h"
00025 #include "radius.h"
00026
00027 void radius_first(void)
00028 {
00029 long int i ,
00030 ip;
00031
00032 bool lgDoPun;
00033
00034 int indexOfSmallest = 0;
00035
00036 const int NUM_DR_TYPES = 13;
00037
00038 struct t_drValues{
00039 double dr;
00040 char whatToSay[40];
00041 } drValues[NUM_DR_TYPES];
00042
00043 double accel,
00044 BigOpacity,
00045 change,
00046 dr912,
00047 drH2 ,
00048 drContPres ,
00049 drOpacity ,
00050 drStromgren,
00051 drTabDen,
00052 dradf,
00053 drcol,
00054 dr_time_dep,
00055 drthrm,
00056 factor,
00057 winddr;
00058 static double drad_last_iteration=-1.;
00059
00060 DEBUG_ENTRY( "radius_first()" );
00061
00062
00063
00064
00065
00066
00067
00068 if( wind.windv > 0. )
00069 {
00070
00071
00072 ASSERT( dense.pden > 0. && dense.wmole > 0. );
00073 accel = 1.3e-10*dense.pden*dense.wmole;
00074 winddr = POW2(wind.windv)/25./accel;
00075 }
00076 else
00077 {
00078 winddr = 1e30;
00079 }
00080
00081
00082 if( StopCalc.taunu > 0.99 && StopCalc.taunu < 3. )
00083 {
00084 dr912 = StopCalc.tauend/6.3e-18/(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac)*Z/50.;
00085 }
00086 else
00087 {
00088 dr912 = 1e30;
00089 }
00090
00091 if( dynamics.lgStatic && iteration > 2 )
00092 {
00093
00094
00095 dr_time_dep = drad_last_iteration;
00096 }
00097 else
00098 {
00099 dr_time_dep = 1e30;
00100 }
00101
00102
00103
00104
00105
00106
00107
00108 if( StopCalc.HColStop < 5e29 )
00109 {
00110
00111 drcol = log10(StopCalc.HColStop) - log10(dense.gas_phase[ipHYDROGEN]*geometry.FillFac* 20.);
00112 }
00113 else if( StopCalc.colpls < 5e29 )
00114 {
00115
00116 drcol = log10(StopCalc.colpls) - log10(dense.xIonDense[ipHYDROGEN][1]*geometry.FillFac* 20.);
00117 }
00118 else if( StopCalc.colnut < 5e29 )
00119 {
00120
00121 drcol = log10(StopCalc.colnut) - log10(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac*50.);
00122 }
00123 else
00124 {
00125
00126 drcol = 30.;
00127 }
00128
00129 drcol = pow(10.,MIN2(35.,drcol));
00130
00131
00132
00133
00134
00135
00136
00137 if( dense.flong != 0. )
00138 {
00139
00140 dradf = 6.283/dense.flong/10.;
00141 dradf = MIN4(dradf,radius.router[iteration-1]*Z,drcol,dr912);
00142 }
00143 else
00144 {
00145 dradf = FLT_MAX;
00146 }
00147
00148
00149
00150
00151 if( (rfield.qhtot>0.) && (rfield.qhtot> rfield.qbal*0.01) && (rfield.uh>1e-10) )
00152 {
00153
00154
00155 drStromgren = (double)(rfield.qhtot)/iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]/
00156 POW2((double)dense.gas_phase[ipHYDROGEN]);
00157
00158
00159 if( drStromgren/radius.rinner > 1. )
00160 {
00161
00162 drStromgren = (double)rfield.qhtot*3./(radius.rinner*
00163 iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]*POW2((double)dense.gas_phase[ipHYDROGEN]) );
00164 drStromgren += 1.;
00165
00166 drStromgren = pow( drStromgren , 0.33333);
00167
00168 drStromgren *= radius.rinner;
00169 }
00170
00171
00172 radius.thickness_stromgren = (realnum)drStromgren;
00173
00174
00175 drStromgren /= 100.;
00176 }
00177 else
00178 {
00179 drStromgren = FLT_MAX;
00180 radius.thickness_stromgren = FLT_MAX;
00181 }
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 ip = ipoint(1e-5);
00193
00194
00195 BigOpacity = 0.;
00196 for( i=ip; i < rfield.nflux; i++ )
00197 {
00198
00199
00200
00201 if( rfield.flux[i]>0. && opac.opacity_abs[i] > BigOpacity )
00202 {
00203 BigOpacity = opac.opacity_abs[i];
00204 }
00205 }
00206
00207
00208
00209
00210
00211 if( BigOpacity > SMALLFLOAT )
00212 {
00213 drOpacity = (radius.drChange/100.)/BigOpacity/geometry.FillFac;
00214 }
00215 else
00216 {
00217 drOpacity = 1e30;
00218 }
00219
00220
00221
00222
00223
00224
00225
00226 drthrm = 1.5e31/MAX2(1.,POW2((double)dense.gas_phase[ipHYDROGEN]));
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00237 {
00238 drTabDen = 1.;
00239 i = 1;
00240 factor = 0.;
00241 while( i < 100 && factor < 0.05 )
00242 {
00243
00244 factor = dense.gas_phase[ipHYDROGEN]/
00245 dense_tabden(radius.Radius+drTabDen, drTabDen );
00246
00247 factor = fabs(factor-1.);
00248 drTabDen *= 2.;
00249 i += 1;
00250 }
00251 drTabDen /= 2.;
00252 }
00253 else
00254 {
00255 drTabDen = 1e30;
00256 }
00257
00258
00259
00260
00261
00262 if( hmi.H2_total/dense.gas_phase[ipHYDROGEN] < 0.1 )
00263 {
00264 change = 0.1;
00265 }
00266 else
00267 {
00268
00269
00270
00271 change = 1.;
00272 }
00273
00274
00275
00276
00277
00278
00279
00280 change = 0.001;
00281
00282 if( h2.lgH2ON && hmi.lgBigH2_evaluated )
00283 {
00284 if( fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > 0.05 )
00285 {
00286
00287
00288
00289
00290 change = 0.0001;
00291 }
00292 else
00293 {
00294
00295
00296 change = 0.001;
00297 }
00298 }
00299 drH2 = change / SDIV(
00300 hmi.H2_total * geometry.FillFac * hmi.H2Opacity );
00301
00302
00303
00304
00305
00306 if( (strcmp( dense.chDenseLaw, "CPRE" )==0) && pressure.lgContRadPresOn )
00307 {
00308
00309 drContPres = 0.05 * pressure.PresTotlCurr /
00310 (wind.AccelTot*dense.xMassDensity*geometry.FillFac*geometry.DirectionalCosin);
00311 }
00312 else if( wind.windv != 0. )
00313 {
00314
00315 double g = fabs(wind.AccelTot-wind.AccelGravity);
00316
00317 drContPres = 0.05*POW2(wind.windv)/(2.*SDIV(g));
00318 }
00319 else
00320 drContPres = 1e30;
00321
00322 drValues[0].dr = drOpacity;
00323 drValues[1].dr = radius.Radius/20.;
00324 drValues[2].dr = drStromgren;
00325 drValues[3].dr = radius.router[iteration-1]/10.;
00326 drValues[4].dr = drcol;
00327 drValues[5].dr = dr912;
00328 drValues[6].dr = drthrm;
00329 drValues[7].dr = winddr;
00330 drValues[8].dr = dradf;
00331 drValues[9].dr = drTabDen;
00332 drValues[10].dr = drH2;
00333 drValues[11].dr = drContPres;
00334 drValues[12].dr = dr_time_dep;
00335
00336 strcpy( drValues[0].whatToSay, "drOpacity" );
00337 strcpy( drValues[1].whatToSay, "radius.Radius/20.");
00338 strcpy( drValues[2].whatToSay, "drStromgren");
00339 strcpy( drValues[3].whatToSay, "radius.router[iteration-1]/10.");
00340 strcpy( drValues[4].whatToSay, "drcol");
00341 strcpy( drValues[5].whatToSay, "dr912");
00342 strcpy( drValues[6].whatToSay, "drthrm");
00343 strcpy( drValues[7].whatToSay, "winddr");
00344 strcpy( drValues[8].whatToSay, "dradf");
00345 strcpy( drValues[9].whatToSay, "drTabDen");
00346 strcpy( drValues[10].whatToSay, "drH2");
00347 strcpy( drValues[11].whatToSay, "drContPres");
00348 strcpy( drValues[12].whatToSay, "dr_time_dep");
00349
00350 for( i=0; i<NUM_DR_TYPES; i++ )
00351 {
00352 if( drValues[i].dr < drValues[indexOfSmallest].dr )
00353 {
00354 indexOfSmallest = i;
00355 }
00356 }
00357
00358 radius.drad = drValues[indexOfSmallest].dr;
00359
00360
00361 if( radius.sdrmin >= radius.drad )
00362 {
00363 radius.drad = radius.sdrmin;
00364
00365
00366
00367 radius.lgDR2Big = true;
00368 }
00369 else
00370 {
00371 radius.lgDR2Big = false;
00372 }
00373
00374
00375
00376 radius.drad = MIN2( radius.sdrmax , radius.drad );
00377
00378 #if 0
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 radius.drad = MIN4( MIN3( drOpacity, radius.Radius/20., drStromgren ),
00389 MIN3( radius.router[iteration-1]/10., drcol, dr912 ),
00390 MIN4( drthrm, winddr, dradf, drTabDen ),
00391 MIN3( drH2, drContPres, dr_time_dep ) );
00392
00393
00394 radius.drad = MAX2( radius.drad , radius.sdrmin );
00395
00396
00397
00398
00399 if( fp_equal( radius.drad, radius.sdrmin ) )
00400 {
00401 radius.lgDR2Big = true;
00402 }
00403 else
00404 {
00405 radius.lgDR2Big = false;
00406 }
00407
00408
00409
00410 radius.drad = MIN2( radius.sdrmax , radius.drad );
00411 #endif
00412
00413 radius.drad_x_fillfac = radius.drad * geometry.FillFac;
00414
00415
00416 drad_last_iteration = radius.drad;
00417
00418
00419
00420 if( radius.lgDrMnOn )
00421 {
00422
00423
00424
00425
00426 radius.drMinimum = (realnum)(radius.drad * dense.gas_phase[ipHYDROGEN]/1e7);
00427 }
00428 else
00429 {
00430 radius.drMinimum = 0.;
00431 }
00432
00433
00434
00435 if( radius.lgSMinON )
00436 {
00437
00438
00439
00440
00441 radius.drMinimum = MIN2(radius.drMinimum * dense.gas_phase[ipHYDROGEN],(realnum)(radius.sdrmin/10.f) );
00442 }
00443
00444 if( trace.lgTrace )
00445 {
00446 fprintf( ioQQQ,
00447 " radius_first called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n",
00448 radius.drad, radius.drMinimum/ dense.gas_phase[ipHYDROGEN], radius.sdrmin, radius.sdrmax );
00449 }
00450
00451 if( radius.drad < SMALLFLOAT*1.1 )
00452 {
00453 fprintf( ioQQQ,
00454 " PROBLEM radius_first detected likely insanity, found dr=%13.5e \n", radius.drad);
00455 fprintf( ioQQQ,
00456 " radius_first: calculation continuing but crash is likely. \n");
00457
00458 TotalInsanity();
00459 }
00460
00461
00462
00463
00464 if( punch.lgDROn )
00465 {
00466 lgDoPun = true;
00467 }
00468 else
00469 {
00470 lgDoPun = false;
00471 }
00472
00473
00474 if( lgDoPun )
00475 {
00476
00477 if( iteration > 1 && punch.lgDRHash )
00478 {
00479 static int iter_punch=-1;
00480 if( iteration !=iter_punch )
00481 fprintf( punch.ipDRout, "%s\n",punch.chHashString );
00482 iter_punch = iteration;
00483 }
00484
00485
00486 fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drad, radius.Depth2Go );
00487
00488 if( radius.lgDR2Big )
00489 {
00490 fprintf( punch.ipDRout,
00491 "radius_first keys from radius.sdrmin\n");
00492
00493 }
00494 else if( fp_equal( radius.drad, radius.sdrmax ) )
00495 {
00496
00497 fprintf( punch.ipDRout,
00498 "radius_first keys from radius.sdrmax\n");
00499 }
00500 else
00501 {
00502 ASSERT( indexOfSmallest < NUM_DR_TYPES - 1 );
00503 fprintf( punch.ipDRout, "radius_first keys from %s\n",
00504 drValues[indexOfSmallest].whatToSay);
00505 }
00506
00507
00508
00509 #if 0
00510 if( fp_equal( radius.drad, drOpacity ) )
00511 {
00512 fprintf( punch.ipDRout,
00513 "radius_first keys from drOpacity, opac was %.2e at %.2e Ryd\n",
00514 BigOpacity , BigOpacityAnu );
00515 }
00516 else if( fp_equal( radius.drad, radius.Radius/20. ) )
00517 {
00518 fprintf( punch.ipDRout,
00519 "radius_first keys from radius.Radius\n" );
00520 }
00521 else if( fp_equal( radius.drad, drStromgren ) )
00522 {
00523 fprintf( punch.ipDRout,
00524 "radius_first keys from drStromgren\n");
00525 }
00526 else if( fp_equal( radius.drad, dr_time_dep ) )
00527 {
00528 fprintf( punch.ipDRout,
00529 "radius_first keys from time dependent\n");
00530 }
00531 else if( fp_equal( radius.drad, radius.router[iteration-1]/10. ) )
00532 {
00533 fprintf( punch.ipDRout,
00534 "radius_first keys from radius.router[iteration-1]\n");
00535 }
00536 else if( fp_equal( radius.drad, drcol ) )
00537 {
00538 fprintf( punch.ipDRout,
00539 "radius_first keys from drcol\n");
00540 }
00541 else if( fp_equal( radius.drad, radius.sdrmin ) )
00542 {
00543 fprintf( punch.ipDRout,
00544 "radius_first keys from radius.sdrmin\n");
00545 }
00546 else if( fp_equal( radius.drad, dr912 ) )
00547 {
00548 fprintf( punch.ipDRout,
00549 "radius_first keys from dr912\n");
00550 }
00551 else if( fp_equal( radius.drad, radius.sdrmax ) )
00552 {
00553 fprintf( punch.ipDRout,
00554 "radius_first keys from radius.sdrmax\n");
00555 }
00556 else if( fp_equal( radius.drad, drthrm ) )
00557 {
00558 fprintf( punch.ipDRout,
00559 "radius_first keys from drthrm\n");
00560 }
00561 else if( fp_equal( radius.drad, winddr ) )
00562 {
00563 fprintf( punch.ipDRout,
00564 "radius_first keys from winddr\n");
00565 }
00566 else if( fp_equal( radius.drad, drH2 ) )
00567 {
00568 fprintf( punch.ipDRout,
00569 "radius_first keys from H2 lyman lines\n");
00570 }
00571 else if( fp_equal( radius.drad, dradf ) )
00572 {
00573 fprintf( punch.ipDRout,
00574 "radius_first keys from dradf\n");
00575 }
00576 else if( fp_equal( radius.drad, drTabDen ) )
00577 {
00578 fprintf( punch.ipDRout,
00579 "radius_first keys from drTabDen\n");
00580 }
00581 else if( fp_equal( radius.drad, drContPres ) )
00582 {
00583 fprintf( punch.ipDRout,
00584 "radius_first keys from radiative acceleration across zone\n");
00585 }
00586 else
00587 {
00588 fprintf( punch.ipDRout, "radius_first insanity\n" );
00589 fprintf( ioQQQ, "radius_first insanity, radius is %e\n" ,
00590 radius.drad);
00591 ShowMe();
00592 }
00593 #endif
00594
00595 }
00596 return;
00597 }