00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "phycon.h"
00009 #include "opacity.h"
00010 #include "rfield.h"
00011 #include "struc.h"
00012 #include "thermal.h"
00013 #include "dense.h"
00014 #include "h2.h"
00015 #include "geometry.h"
00016 #include "conv.h"
00017 #include "dynamics.h"
00018 #include "radius.h"
00019 #include "zones.h"
00020 #include "doppvel.h"
00021
00022 #define IOFF 3
00023
00024 void ZoneStart(const char *chMode)
00025 {
00026 bool lgNeOscil,
00027 lgTeOscil;
00028 long int i;
00029
00030 double dt1 , de1,
00031
00032
00033
00034 DilutionCorrec ,
00035 drFac ,
00036 dTauThisZone,
00037 outwrd,
00038 ratio,
00039 rm,
00040 rad_middle_zone,
00041 vin,
00042 vout;
00043
00044 static double DTaver , DEaver,
00045 dt2 , de2;
00046
00047 DEBUG_ENTRY( "ZoneStart()" );
00048
00050
00051 if( DoppVel.lgTurbLawOn )
00052 {
00053 DoppVel.TurbVel = DoppVel.TurbVelZero *
00054 pow( (realnum)(radius.Radius/radius.rinner), DoppVel.TurbVelLaw );
00055 }
00056
00057
00058
00059
00060
00061
00062
00063
00064 geometry.FillFac = (realnum)(geometry.fiscal*
00065 pow( radius.Radius/radius.rinner ,(double)geometry.filpow));
00066 geometry.FillFac = (realnum)MIN2(1.,geometry.FillFac);
00067
00068 if( strcmp(chMode,"init") == 0 )
00069 {
00070
00071
00072
00073
00074 radius.depth = 1e-30;
00075
00076
00077 radius.depth_x_fillfac = radius.depth*geometry.FillFac;
00078
00079 radius.drad_x_fillfac = radius.drad*geometry.FillFac;
00080
00081
00082 radius.Radius = radius.rinner;
00083 radius.Radius_mid_zone = radius.rinner + radius.dRadSign*radius.drad/2.;
00084
00085
00086 radius.drNext = radius.drad;
00087
00088
00089 radius.depth_mid_zone = radius.drad/2.;
00090
00091
00092 radius.glbdst = radius.glbrad;
00093
00094
00095 if( radius.StopThickness[iteration-1] < 5e29 )
00096 {
00097 radius.Depth2Go = radius.StopThickness[iteration-1];
00098 }
00099 else if( iteration > 1 )
00100 {
00101
00102 radius.Depth2Go = radius.StopThickness[iteration-2];
00103 }
00104 else
00105 {
00106
00107 radius.Depth2Go = 0.;
00108 }
00109 }
00110 else if( strcmp(chMode,"incr") == 0 )
00111 {
00112
00113 radius.drad_mid_zone = (radius.drad+radius.drNext)/2.;
00114 radius.drad = radius.drNext;
00115
00116
00117 radius.dr_min_last_iter = MIN2( radius.dr_min_last_iter , radius.drNext );
00118 radius.dr_max_last_iter = MAX2( radius.dr_max_last_iter , radius.drNext );
00119
00120 ASSERT( radius.drad > 0. );
00121 radius.depth += radius.drad;
00122 radius.depth_mid_zone = radius.depth - radius.drad/2.;
00123
00124
00125 radius.drad_x_fillfac = radius.drad*geometry.FillFac;
00126
00127
00128 radius.depth_x_fillfac += radius.drad_x_fillfac;
00129
00130
00131 radius.Depth2Go = MAX2( 0.,radius.Depth2Go - radius.drad );
00132
00133
00134
00135 radius.Radius += radius.drad*radius.dRadSign;
00136
00137
00138
00139 radius.Radius_mid_zone = radius.Radius - radius.dRadSign*radius.drad/2.;
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 drFac = radius.drad*geometry.FillFac*geometry.DirectionalCosin*1.35;
00163
00164
00165
00166
00167 DilutionCorrec = 1./POW2(
00168 (radius.Radius-radius.dRadSign*radius.drad/2.)/(radius.Radius-radius.dRadSign*radius.drad) );
00169
00170
00171
00172
00173
00174
00175 for( i=0; i <= rfield.nflux; i++ )
00176 {
00177 dTauThisZone = opac.opacity_abs[i]*drFac;
00178
00179
00180
00181
00182
00183 if( dTauThisZone < 1e-4 )
00184 {
00185
00186 opac.tmn[i] = 1.f;
00187 }
00188 else if( dTauThisZone < 5. )
00189 {
00190
00191 opac.tmn[i] = (realnum)((1. - exp(-dTauThisZone))/(dTauThisZone));
00192 }
00193 else
00194 {
00195
00196 opac.tmn[i] = (realnum)(1./(dTauThisZone));
00197 }
00198
00199
00200 opac.tmn[i] *= (realnum)DilutionCorrec;
00201
00202 rfield.flux_beam_const[i] *= opac.tmn[i];
00203 rfield.flux_beam_time[i] *= opac.tmn[i];
00204 rfield.flux_isotropic[i] *= opac.tmn[i];
00205 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
00206 rfield.flux_isotropic[i];
00207
00208
00209 rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
00210 }
00211
00212
00213 radius.glbdst -= (realnum)radius.drad;
00214
00215
00216
00217
00218
00219
00220 if( (strcmp(dense.chDenseLaw,"CPRE") != 0) &&
00221 thermal.lgPredNextTe &&
00222 (nzone > IOFF + 1) )
00223 {
00224 phycon.TeInit = phycon.te;
00225 phycon.EdenInit = dense.eden;
00226 lgTeOscil = false;
00227 lgNeOscil = false;
00228 dt1 = 0.;
00229 dt2 = 0.;
00230 de1 = 0.;
00231 de2 = 0.;
00232 DTaver = 0.;
00233 DEaver = 0.;
00234 for( i=nzone - IOFF-1; i < (nzone - 1); i++ )
00235 {
00236 dt1 = dt2;
00237 de1 = de2;
00238
00239
00240 dt2 = struc.testr[i-1] - struc.testr[i];
00241 de2 = struc.ednstr[i-1] - struc.ednstr[i];
00242 DTaver += dt2;
00243 DEaver += de2;
00244 if( dt1*dt2 < 0. )
00245 {
00246 lgTeOscil = true;
00247 }
00248 if( de1*de2 < 0. )
00249 {
00250 lgNeOscil = true;
00251 }
00252 }
00253
00254
00255 if( !lgTeOscil )
00256 {
00257 DTaver /= (double)(IOFF);
00258
00259 dt2 = fabs(dt2);
00260 rm = fabs(DTaver);
00261 DTaver = sign(MIN2(rm,dt2),DTaver);
00262
00263
00264
00265
00266 DTaver = sign(MIN2(rm,0.01*phycon.te),DTaver);
00267
00268 TempChange(phycon.te - DTaver , true);
00269 }
00270 else
00271 {
00272
00273 DTaver = 0.;
00274 }
00275
00276
00277 phycon.TeProp = phycon.te;
00278
00279
00280 if( !lgNeOscil )
00281 {
00282 DEaver /= IOFF;
00283 de2 = fabs(de2);
00284 rm = fabs(DEaver);
00285
00286 DEaver = sign(MIN2(rm,de2),DEaver);
00287
00288 DEaver = sign(MIN2(rm,0.05*dense.eden),DEaver);
00289
00290 dense.eden -= DEaver;
00291 }
00292 else
00293 {
00294
00295 DEaver = 0.;
00296 }
00297
00298
00299 phycon.EdenProp = dense.eden;
00300
00301
00302 TempChange(phycon.te , false);
00303 }
00304 }
00305
00306 else
00307 {
00308 fprintf( ioQQQ, " PROBLEM ZoneStart called with insane argument, %4.4s\n",
00309 chMode );
00310
00311 TotalInsanity();
00312 }
00313
00314
00315
00316
00317 radius.drad_x_fillfac_mean = radius.drad_x_fillfac;
00318 i = 1;
00319 while( i < 5 && nzone-i-1>=0 )
00320 {
00321 radius.drad_x_fillfac_mean += struc.drad_x_fillfac[nzone-i-1];
00322 ++i;
00323 }
00324 radius.drad_x_fillfac_mean /= i;
00325
00326
00327 if( dynamics.lgAdvection )
00328 DynaStartZone();
00329
00330
00331
00332
00333
00334
00335
00336 h2.nCallH2_this_zone = 0;
00337
00338
00339 if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00340 {
00341 ratio = radius.drad/(radius.glbdst + radius.drad);
00342
00343 if( ratio < 1e-14 )
00344 {
00345 radius.lgdR2Small = true;
00346 }
00347 else
00348 {
00349 radius.lgdR2Small = false;
00350 }
00351 }
00352
00353
00354
00355
00356
00357
00358
00359 radius.r1r0sq = POW2( radius.Radius/radius.rinner );
00360
00361
00362 radius.dRNeff = radius.drNext*geometry.FillFac;
00363
00364
00365 rad_middle_zone = radius.Radius - radius.dRadSign*radius.drad/2.;
00366
00367
00368 if( radius.drad/radius.Radius > 0.01 )
00369 {
00370 double Ropp = radius.Radius - radius.dRadSign*radius.drad;
00371 radius.darea_x_fillfac = PI*abs(pow2(radius.Radius) - pow2(Ropp))*geometry.FillFac;
00372 }
00373 else
00374 radius.darea_x_fillfac = PI2*rad_middle_zone*radius.drad*geometry.FillFac;
00375
00376
00377
00378
00379
00380 if( radius.drad/radius.Radius > 0.01 )
00381 {
00382 double r1 = radius.Radius - radius.dRadSign*radius.drad;
00383 double rin_zone = min(r1,radius.Radius);
00384 double rout_zone = max(r1,radius.Radius);
00385 vin = pow2(rin_zone/radius.rinner)*rin_zone/3.;
00386 if( rin_zone > radius.CylindHigh )
00387 {
00388
00389
00390
00391 double h = rin_zone-radius.CylindHigh;
00392 double v2cap = pow2(h/radius.rinner)*(rin_zone - h/3.)/2.;
00393 vin -= v2cap;
00394 }
00395 vout = pow2(rout_zone/radius.rinner)*rout_zone/3.;
00396 if( rout_zone > radius.CylindHigh )
00397 {
00398 double h = rout_zone-radius.CylindHigh;
00399 double v2cap = pow2(h/radius.rinner)*(rout_zone - h/3.)/2.;
00400 vout -= v2cap;
00401 }
00402
00403
00404 radius.dVeffVol = (vout - vin)*geometry.FillFac;
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415 if( geometry.iEmissPower == 2 )
00416 {
00417 radius.dVeffAper = radius.dVeffVol;
00418 }
00419 else if( geometry.iEmissPower == 1 )
00420 {
00421 double ain = (rin_zone/radius.rinner)*rin_zone/2.;
00422 if( rin_zone > radius.CylindHigh )
00423 {
00424
00425
00426
00427 double Theta = 2.*acos(min(radius.CylindHigh/rin_zone,1.));
00428 ain *= 1. - max(Theta - sin(Theta),0.)/PI;
00429 }
00430 double aout = (rout_zone/radius.rinner)*rout_zone/2.;
00431 if( rout_zone > radius.CylindHigh )
00432 {
00433 double Theta = 2.*acos(min(radius.CylindHigh/rout_zone,1.));
00434 aout *= 1. - max(Theta - sin(Theta),0.)/PI;
00435 }
00436 radius.dVeffAper = (aout - ain)*geometry.FillFac;
00437 }
00438 else if( geometry.iEmissPower == 0 )
00439 {
00440 radius.dVeffAper = radius.drad*geometry.FillFac;
00441 }
00442 }
00443 else
00444 {
00445
00446
00447 radius.dVeffVol = (rad_middle_zone/radius.rinner)*radius.drad*geometry.FillFac*
00448 (min(rad_middle_zone,radius.CylindHigh)/radius.rinner);
00449 if( geometry.iEmissPower == 2 )
00450 {
00451 radius.dVeffAper = radius.dVeffVol;
00452 }
00453 else if( geometry.iEmissPower == 1 )
00454 {
00455 radius.dVeffAper = (rad_middle_zone/radius.rinner)*radius.drad*geometry.FillFac;
00456 if( rad_middle_zone > radius.CylindHigh )
00457 {
00458 double Theta = 2.*acos(min(radius.CylindHigh/rad_middle_zone,1.));
00459 double q = sqrt(max(1.-pow2(radius.CylindHigh/rad_middle_zone),0.))*rad_middle_zone;
00460 radius.dVeffAper *= 1. - max(Theta - sin(Theta),0.)/PI - max(1. - cos(Theta),0.)*radius.CylindHigh/(PI*q);
00461 }
00462 }
00463 else if( geometry.iEmissPower == 0 )
00464 {
00465 radius.dVeffAper = radius.drad*geometry.FillFac;
00466 }
00467 }
00468
00469
00470 radius.dVeffVol *= geometry.covgeo;
00471 radius.dVeffAper *= geometry.covaper;
00472
00473
00474
00475
00476 outwrd = (1. + geometry.covrt)/2.;
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 radius.dVolOutwrd = outwrd*radius.drad_x_fillfac;
00488 ASSERT( radius.dVolOutwrd > 0. );
00489
00490
00491
00492
00493
00494
00495
00496 radius.dVolReflec = (1. - outwrd)*radius.drad_x_fillfac*radius.r1r0sq;
00497
00498 if( geometry.lgSphere )
00499 {
00500
00501 radius.BeamInIn = 0.;
00502 radius.BeamInOut = radius.Radius*radius.drad_x_fillfac/(radius.Radius +
00503 2.*radius.drad);
00504 }
00505 else
00506 {
00507 radius.BeamInIn = radius.drad_x_fillfac*radius.r1r0sq;
00508
00509
00510 radius.BeamInOut = 0.;
00511 }
00512
00513 radius.BeamOutOut = radius.Radius*radius.drad_x_fillfac/(radius.Radius +
00514 2.*radius.drad);
00515 return;
00516 }
00517
00518 void ZoneEnd(void)
00519 {
00520 long i;
00521
00522 DEBUG_ENTRY( "ZoneEnd()" );
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536 for( i=0; i <= rfield.nflux; i++ )
00537 {
00538 rfield.flux_beam_const[i] /= opac.tmn[i];
00539 rfield.flux_beam_time[i] /= opac.tmn[i];
00540 rfield.flux_isotropic[i] /= opac.tmn[i];
00541 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
00542 rfield.flux_isotropic[i];
00543
00544 rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
00545 }
00546
00547
00548 if( dynamics.lgAdvection )
00549 DynaEndZone();
00550 return;
00551 }