33 int indexOfSmallest = 0;
35 const double Z = 1.0001;
36 const int NUM_DR_TYPES = 13;
41 } drValues[NUM_DR_TYPES];
58 static double drad_last_iteration=-1.;
68 if(
wind.lgBallistic() )
95 dr_time_dep = drad_last_iteration;
129 drcol = pow(10.,
MIN2(35.,drcol));
137 if(
dense.flong != 0. )
140 dradf = 6.283/
dense.flong/10.;
159 drStromgren =
MIN2(1e28 , drStromgren );
162 if( drStromgren/
radius.rinner > 1. )
169 drStromgren = pow( drStromgren , 0.33333);
171 drStromgren *=
radius.rinner;
183 drStromgren = FLT_MAX;
184 radius.thickness_stromgren = FLT_MAX;
200 for( i=ip; i <
rfield.nflux; i++ )
205 if(
rfield.flux[0][i]>0. &&
opac.opacity_abs[i] > BigOpacity )
207 BigOpacity =
opac.opacity_abs[i];
233 drthrm =
MAX2( 0.15, drthrm );
243 if( strcmp(
dense.chDenseLaw,
"DLW2") == 0 )
248 while( i < 100 && factor < 0.05 &&
radius.Radius+drTabDen*2.<
radius.StopThickness[0] )
254 factor = fabs(factor-1.);
289 if(
h2.lgEnabled &&
h2.lgEvaluated )
291 if( fabs(
h2.HeatDexc)/
thermal.ctot > 0.05 )
306 drH2 = change /
SDIV(
313 if( (strcmp(
dense.chDenseLaw,
"CPRE" )==0) &&
pressure.lgContRadPresOn )
316 drContPres = 0.05 *
pressure.PresTotlCurr /
319 else if( !
wind.lgStatic() )
322 double g = fabs(
wind.AccelTotalOutward-
wind.AccelGravity);
329 drValues[0].dr = drOpacity;
330 drValues[1].dr =
radius.Radius/20.;
331 drValues[2].dr = drStromgren;
333 drValues[4].dr = drcol;
334 drValues[5].dr = dr912;
335 drValues[6].dr = drthrm;
336 drValues[7].dr = winddr;
337 drValues[8].dr = dradf;
338 drValues[9].dr = drTabDen;
339 drValues[10].dr = drH2;
340 drValues[11].dr = drContPres;
341 drValues[12].dr = dr_time_dep;
343 strcpy( drValues[0].whatToSay,
"drOpacity" );
344 strcpy( drValues[1].whatToSay,
"radius.Radius/20.");
345 strcpy( drValues[2].whatToSay,
"drStromgren");
346 strcpy( drValues[3].whatToSay,
"radius.StopThickness[iteration-1]/10.");
347 strcpy( drValues[4].whatToSay,
"drcol");
348 strcpy( drValues[5].whatToSay,
"dr912");
349 strcpy( drValues[6].whatToSay,
"drthrm");
350 strcpy( drValues[7].whatToSay,
"winddr");
351 strcpy( drValues[8].whatToSay,
"dradf");
352 strcpy( drValues[9].whatToSay,
"drTabDen");
353 strcpy( drValues[10].whatToSay,
"drH2");
354 strcpy( drValues[11].whatToSay,
"drContPres");
355 strcpy( drValues[12].whatToSay,
"dr_time_dep");
357 for( i=0; i<NUM_DR_TYPES; i++ )
359 if( drValues[i].dr < drValues[indexOfSmallest].dr )
365 radius.drad = drValues[indexOfSmallest].dr;
367 double rfacmin =
radius.lgSdrminRel ?
radius.Radius : 1.;
385 double rfacmax =
radius.lgSdrmaxRel ?
radius.Radius : 1.;
401 MIN4( drthrm, winddr, dradf, drTabDen ),
402 MIN3( drH2, drContPres, dr_time_dep ) );
427 drad_last_iteration =
radius.drad;
459 " radius_first called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n",
467 " PROBLEM radius_first detected likely insanity, found dr=%13.5e \n",
radius.drad);
469 " radius_first: calculation continuing but crash is likely. \n");
492 static int iter_punch=-1;
494 fprintf(
save.ipDRout,
"%s\n",
save.chHashString );
503 fprintf(
save.ipDRout,
504 "radius_first keys from radius.sdrmin\n");
510 fprintf(
save.ipDRout,
511 "radius_first keys from radius.sdrmax\n");
515 ASSERT( indexOfSmallest < NUM_DR_TYPES - 1 );
516 fprintf(
save.ipDRout,
"radius_first keys from %s\n",
517 drValues[indexOfSmallest].whatToSay);
525 fprintf(
save.ipDRout,
526 "radius_first keys from drOpacity, opac was %.2e at %.2e Ryd\n",
527 BigOpacity , BigOpacityAnu );
531 fprintf(
save.ipDRout,
532 "radius_first keys from radius.Radius\n" );
536 fprintf(
save.ipDRout,
537 "radius_first keys from drStromgren\n");
541 fprintf(
save.ipDRout,
542 "radius_first keys from time dependent\n");
546 fprintf(
save.ipDRout,
547 "radius_first keys from radius.StopThickness[iteration-1]\n");
551 fprintf(
save.ipDRout,
552 "radius_first keys from drcol\n");
556 fprintf(
save.ipDRout,
557 "radius_first keys from radius.sdrmin\n");
561 fprintf(
save.ipDRout,
562 "radius_first keys from dr912\n");
566 fprintf(
save.ipDRout,
567 "radius_first keys from radius.sdrmax\n");
571 fprintf(
save.ipDRout,
572 "radius_first keys from drthrm\n");
576 fprintf(
save.ipDRout,
577 "radius_first keys from winddr\n");
581 fprintf(
save.ipDRout,
582 "radius_first keys from H2 lyman lines\n");
586 fprintf(
save.ipDRout,
587 "radius_first keys from dradf\n");
591 fprintf(
save.ipDRout,
592 "radius_first keys from drTabDen\n");
596 fprintf(
save.ipDRout,
597 "radius_first keys from radiative acceleration across zone\n");
601 fprintf(
save.ipDRout,
"radius_first insanity\n" );
602 fprintf(
ioQQQ,
"radius_first insanity, radius is %e\n" ,