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