00001
00002
00003
00004
00005
00006
00007
00008 #include "cddefines.h"
00009 #include "wind.h"
00010 #include "stopcalc.h"
00011 #include "thermal.h"
00012 #include "dynamics.h"
00013 #include "trace.h"
00014 #include "save.h"
00015 #include "pressure.h"
00016 #include "iso.h"
00017 #include "h2.h"
00018 #include "rfield.h"
00019 #include "dense.h"
00020 #include "hmi.h"
00021 #include "geometry.h"
00022 #include "opacity.h"
00023 #include "ipoint.h"
00024 #include "radius.h"
00025
00026 void radius_first(void)
00027 {
00028 long int i ,
00029 ip;
00030
00031 bool lgDoPun;
00032
00033 int indexOfSmallest = 0;
00034
00035 const double Z = 1.0001;
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.lgBallistic() )
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.lgTimeDependentStatic && 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.StopThickness[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_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB/
00156 POW2((double)dense.gas_phase[ipHYDROGEN]);
00157
00158
00159 drStromgren = MIN2(1e28 , drStromgren );
00160
00161
00162 if( drStromgren/radius.rinner > 1. )
00163 {
00164
00165 drStromgren = (double)rfield.qhtot*3./(radius.rinner*
00166 iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB*POW2((double)dense.gas_phase[ipHYDROGEN]) );
00167 drStromgren += 1.;
00168
00169 drStromgren = pow( drStromgren , 0.33333);
00170
00171 drStromgren *= radius.rinner;
00172 }
00173
00174
00175 radius.thickness_stromgren = (realnum)drStromgren;
00176
00177
00178 if( drStromgren > SMALLFLOAT *100.)
00179 drStromgren /= 100.;
00180 }
00181 else
00182 {
00183 drStromgren = FLT_MAX;
00184 radius.thickness_stromgren = FLT_MAX;
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 ip = ipoint(1e-5);
00197
00198
00199 BigOpacity = 0.;
00200 for( i=ip; i < rfield.nflux; i++ )
00201 {
00202
00203
00204
00205 if( rfield.flux[0][i]>0. && opac.opacity_abs[i] > BigOpacity )
00206 {
00207 BigOpacity = opac.opacity_abs[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 drthrm = MAX2( 0.15, drthrm );
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00244 {
00245 drTabDen = 1.;
00246 i = 1;
00247 factor = 0.;
00248 while( i < 100 && factor < 0.05 && radius.Radius+drTabDen*2.<radius.StopThickness[0] )
00249 {
00250
00251 factor = dense.gas_phase[ipHYDROGEN]/
00252 dense_tabden(radius.Radius+drTabDen, drTabDen );
00253
00254 factor = fabs(factor-1.);
00255 drTabDen *= 2.;
00256 i += 1;
00257 }
00258 drTabDen /= 2.;
00259 }
00260 else
00261 {
00262 drTabDen = 1e30;
00263 }
00264
00265
00266
00267
00268
00269 if( hmi.H2_total/dense.gas_phase[ipHYDROGEN] < 0.1 )
00270 {
00271 change = 0.1;
00272 }
00273 else
00274 {
00275
00276
00277
00278 change = 1.;
00279 }
00280
00281
00282
00283
00284
00285
00286
00287 change = 0.001;
00288
00289 if( h2.lgEnabled && h2.lgEvaluated )
00290 {
00291 if( fabs(h2.HeatDexc)/thermal.ctot > 0.05 )
00292 {
00293
00294
00295
00296
00297 change = 0.0001;
00298 }
00299 else
00300 {
00301
00302
00303 change = 0.001;
00304 }
00305 }
00306 drH2 = change / SDIV(
00307 hmi.H2_total * geometry.FillFac * hmi.H2Opacity );
00308
00309
00310
00311
00312
00313 if( (strcmp( dense.chDenseLaw, "CPRE" )==0) && pressure.lgContRadPresOn )
00314 {
00315
00316 drContPres = 0.05 * pressure.PresTotlCurr /
00317 ((double)wind.AccelTotalOutward*dense.xMassDensity*geometry.FillFac*geometry.DirectionalCosin);
00318 }
00319 else if( !wind.lgStatic() )
00320 {
00321
00322 double g = fabs(wind.AccelTotalOutward-wind.AccelGravity);
00323
00324 drContPres = 0.05*POW2(wind.windv)/(2.*SDIV(g));
00325 }
00326 else
00327 drContPres = 1e30;
00328
00329 drValues[0].dr = drOpacity;
00330 drValues[1].dr = radius.Radius/20.;
00331 drValues[2].dr = drStromgren;
00332 drValues[3].dr = radius.StopThickness[iteration-1]/10.;
00333 drValues[4].dr = drcol;
00334 drValues[5].dr = dr912;
00335 drValues[6].dr = drthrm;
00336 drValues[7].dr = winddr;
00337 drValues[8].dr = dradf;
00338 drValues[9].dr = drTabDen;
00339 drValues[10].dr = drH2;
00340 drValues[11].dr = drContPres;
00341 drValues[12].dr = dr_time_dep;
00342
00343 strcpy( drValues[0].whatToSay, "drOpacity" );
00344 strcpy( drValues[1].whatToSay, "radius.Radius/20.");
00345 strcpy( drValues[2].whatToSay, "drStromgren");
00346 strcpy( drValues[3].whatToSay, "radius.StopThickness[iteration-1]/10.");
00347 strcpy( drValues[4].whatToSay, "drcol");
00348 strcpy( drValues[5].whatToSay, "dr912");
00349 strcpy( drValues[6].whatToSay, "drthrm");
00350 strcpy( drValues[7].whatToSay, "winddr");
00351 strcpy( drValues[8].whatToSay, "dradf");
00352 strcpy( drValues[9].whatToSay, "drTabDen");
00353 strcpy( drValues[10].whatToSay, "drH2");
00354 strcpy( drValues[11].whatToSay, "drContPres");
00355 strcpy( drValues[12].whatToSay, "dr_time_dep");
00356
00357 for( i=0; i<NUM_DR_TYPES; i++ )
00358 {
00359 if( drValues[i].dr < drValues[indexOfSmallest].dr )
00360 {
00361 indexOfSmallest = i;
00362 }
00363 }
00364
00365 radius.drad = drValues[indexOfSmallest].dr;
00366
00367 double rfacmin = radius.lgSdrminRel ? radius.Radius : 1.;
00368
00369 if( rfacmin*radius.sdrmin >= radius.drad )
00370 {
00371 radius.drad = rfacmin*radius.sdrmin;
00372
00373
00374
00375 radius.lgDR2Big = true;
00376 }
00377 else
00378 {
00379 radius.lgDR2Big = false;
00380 }
00381
00382
00383
00384
00385 double rfacmax = radius.lgSdrmaxRel ? radius.Radius : 1.;
00386 radius.drad = MIN2( rfacmax*radius.sdrmax, radius.drad );
00387 radius.drad_mid_zone = radius.drad/2.;
00388
00389 #if 0
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399 radius.drad = MIN4( MIN3( drOpacity, radius.Radius/20., drStromgren ),
00400 MIN3( radius.StopThickness[iteration-1]/10., drcol, dr912 ),
00401 MIN4( drthrm, winddr, dradf, drTabDen ),
00402 MIN3( drH2, drContPres, dr_time_dep ) );
00403
00404
00405 radius.drad = MAX2( radius.drad, rfacmin*radius.sdrmin );
00406
00407
00408
00409
00410 if( fp_equal( radius.drad, rfacmin*radius.sdrmin ) )
00411 {
00412 radius.lgDR2Big = true;
00413 }
00414 else
00415 {
00416 radius.lgDR2Big = false;
00417 }
00418
00419
00420
00421 radius.drad = MIN2( rfacmax*radius.sdrmax, radius.drad );
00422 #endif
00423
00424 radius.drad_x_fillfac = radius.drad * geometry.FillFac;
00425
00426
00427 drad_last_iteration = radius.drad;
00428
00429
00430
00431 if( radius.lgDrMnOn )
00432 {
00433
00434
00435
00436
00437 radius.drMinimum = (realnum)(radius.drad * dense.gas_phase[ipHYDROGEN]/1e7);
00438 }
00439 else
00440 {
00441 radius.drMinimum = 0.;
00442 }
00443
00444
00445
00446 if( radius.lgSMinON )
00447 {
00448
00449
00450
00451
00452 radius.drMinimum = MIN2(radius.drMinimum * dense.gas_phase[ipHYDROGEN],
00453 (realnum)(rfacmin*radius.sdrmin/10.f) );
00454 }
00455
00456 if( trace.lgTrace )
00457 {
00458 fprintf( ioQQQ,
00459 " radius_first called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n",
00460 radius.drad, radius.drMinimum/ dense.gas_phase[ipHYDROGEN],
00461 rfacmin*radius.sdrmin, rfacmax*radius.sdrmax );
00462 }
00463
00464 if( radius.drad < SMALLFLOAT*1.1 )
00465 {
00466 fprintf( ioQQQ,
00467 " PROBLEM radius_first detected likely insanity, found dr=%13.5e \n", radius.drad);
00468 fprintf( ioQQQ,
00469 " radius_first: calculation continuing but crash is likely. \n");
00470
00471 TotalInsanity();
00472 }
00473
00474
00475
00476
00477 if( save.lgDROn )
00478 {
00479 lgDoPun = true;
00480 }
00481 else
00482 {
00483 lgDoPun = false;
00484 }
00485
00486
00487 if( lgDoPun )
00488 {
00489
00490 if( iteration > 1 && save.lgDRHash )
00491 {
00492 static int iter_punch=-1;
00493 if( iteration !=iter_punch )
00494 fprintf( save.ipDRout, "%s\n",save.chHashString );
00495 iter_punch = iteration;
00496 }
00497
00498
00499 fprintf( save.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drad, radius.Depth2Go );
00500
00501 if( radius.lgDR2Big )
00502 {
00503 fprintf( save.ipDRout,
00504 "radius_first keys from radius.sdrmin\n");
00505
00506 }
00507 else if( fp_equal( radius.drad, rfacmax*radius.sdrmax ) )
00508 {
00509
00510 fprintf( save.ipDRout,
00511 "radius_first keys from radius.sdrmax\n");
00512 }
00513 else
00514 {
00515 ASSERT( indexOfSmallest < NUM_DR_TYPES - 1 );
00516 fprintf( save.ipDRout, "radius_first keys from %s\n",
00517 drValues[indexOfSmallest].whatToSay);
00518 }
00519
00520
00521
00522 #if 0
00523 if( fp_equal( radius.drad, drOpacity ) )
00524 {
00525 fprintf( save.ipDRout,
00526 "radius_first keys from drOpacity, opac was %.2e at %.2e Ryd\n",
00527 BigOpacity , BigOpacityAnu );
00528 }
00529 else if( fp_equal( radius.drad, radius.Radius/20. ) )
00530 {
00531 fprintf( save.ipDRout,
00532 "radius_first keys from radius.Radius\n" );
00533 }
00534 else if( fp_equal( radius.drad, drStromgren ) )
00535 {
00536 fprintf( save.ipDRout,
00537 "radius_first keys from drStromgren\n");
00538 }
00539 else if( fp_equal( radius.drad, dr_time_dep ) )
00540 {
00541 fprintf( save.ipDRout,
00542 "radius_first keys from time dependent\n");
00543 }
00544 else if( fp_equal( radius.drad, radius.StopThickness[iteration-1]/10. ) )
00545 {
00546 fprintf( save.ipDRout,
00547 "radius_first keys from radius.StopThickness[iteration-1]\n");
00548 }
00549 else if( fp_equal( radius.drad, drcol ) )
00550 {
00551 fprintf( save.ipDRout,
00552 "radius_first keys from drcol\n");
00553 }
00554 else if( fp_equal( radius.drad, rfacmin*radius.sdrmin ) )
00555 {
00556 fprintf( save.ipDRout,
00557 "radius_first keys from radius.sdrmin\n");
00558 }
00559 else if( fp_equal( radius.drad, dr912 ) )
00560 {
00561 fprintf( save.ipDRout,
00562 "radius_first keys from dr912\n");
00563 }
00564 else if( fp_equal( radius.drad, rfacmax*radius.sdrmax ) )
00565 {
00566 fprintf( save.ipDRout,
00567 "radius_first keys from radius.sdrmax\n");
00568 }
00569 else if( fp_equal( radius.drad, drthrm ) )
00570 {
00571 fprintf( save.ipDRout,
00572 "radius_first keys from drthrm\n");
00573 }
00574 else if( fp_equal( radius.drad, winddr ) )
00575 {
00576 fprintf( save.ipDRout,
00577 "radius_first keys from winddr\n");
00578 }
00579 else if( fp_equal( radius.drad, drH2 ) )
00580 {
00581 fprintf( save.ipDRout,
00582 "radius_first keys from H2 lyman lines\n");
00583 }
00584 else if( fp_equal( radius.drad, dradf ) )
00585 {
00586 fprintf( save.ipDRout,
00587 "radius_first keys from dradf\n");
00588 }
00589 else if( fp_equal( radius.drad, drTabDen ) )
00590 {
00591 fprintf( save.ipDRout,
00592 "radius_first keys from drTabDen\n");
00593 }
00594 else if( fp_equal( radius.drad, drContPres ) )
00595 {
00596 fprintf( save.ipDRout,
00597 "radius_first keys from radiative acceleration across zone\n");
00598 }
00599 else
00600 {
00601 fprintf( save.ipDRout, "radius_first insanity\n" );
00602 fprintf( ioQQQ, "radius_first insanity, radius is %e\n" ,
00603 radius.drad);
00604 ShowMe();
00605 }
00606 #endif
00607
00608 }
00609 return;
00610 }