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