00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "cddefines.h"
00011 #include "abund.h"
00012 #include "hmi.h"
00013 #include "struc.h"
00014 #include "trace.h"
00015 #include "wind.h"
00016 #include "phycon.h"
00017 #include "thermal.h"
00018 #include "dense.h"
00019 #include "geometry.h"
00020 #include "radius.h"
00021 #include "mole.h"
00022 #include "dynamics.h"
00023 #include "pressure.h"
00024 #include "colden.h"
00025 #include "conv.h"
00026 #include "cosmology.h"
00027 #include "dark_matter.h"
00028
00029
00030 static double pressure_change_factor;
00031
00032
00033
00034
00035 int PressureChange(
00036
00037 double dP_chng_factor )
00038 {
00039 long int ion,
00040 nelem,
00041 lgChange,
00042 mol;
00043
00044 double abun,
00045 edensave,
00046 hold;
00047
00048
00049
00050 double pdelta;
00051
00052 static double FacAbun,
00053 FacAbunSav,
00054 OldAbun;
00055
00056 DEBUG_ENTRY( "PressureChange()" );
00057
00058 edensave = dense.eden;
00059
00060
00061
00062
00063 PresTotCurrent();
00064
00065
00066
00067 if( nzone != conv.hist_pres_nzone )
00068 {
00069
00070 conv.hist_pres_nzone = nzone;
00071 conv.hist_pres_density.clear();
00072 conv.hist_pres_current.clear();
00073 conv.hist_pres_correct.clear();
00074 }
00075
00076
00077 conv.hist_pres_density.push_back( dense.gas_phase[ipHYDROGEN] );
00078 conv.hist_pres_current.push_back( pressure.PresTotlCurr );
00079 conv.hist_pres_correct.push_back( pressure.PresTotlCorrect );
00080
00081
00082 hold = dense.gas_phase[ipHYDROGEN];
00083
00084
00085 lgChange = false;
00086
00087
00088
00089
00090
00091 conv.lgConvPres = lgConvPres();
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 if( pressure_change_factor != 1. )
00107 {
00108 lgChange = true;
00109 }
00110
00111
00112
00113 pdelta = 0.03 * dP_chng_factor;
00114
00115
00116 pressure_change_factor = MIN2(pressure_change_factor,1.+pdelta);
00117 pressure_change_factor = MAX2(pressure_change_factor,1.-pdelta);
00118
00119 {
00120
00121 enum{DEBUG_LOC=false};
00122 static long int nsave=-1;
00123
00124 if( DEBUG_LOC )
00125 {
00126 if( nsave-nzone ) fprintf(ioQQQ,"\n");
00127 nsave = nzone;
00128 fprintf(ioQQQ,"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\n",
00129 nzone,
00130 pressure_change_factor,
00131
00132 (pressure.PresTotlCorrect-pressure.PresTotlCurr)*100./pressure.PresTotlCorrect,
00133 pressure.PresTotlCorrect,
00134 pressure.PresTotlCurr,
00135 pressure.PresGasCurr,
00136 pressure.pres_radiation_lines_curr );
00137 }
00138 }
00139
00140
00141
00142 if( abund.lgAbTaON )
00143 {
00144 lgChange = true;
00145 for( nelem=1; nelem < LIMELM; nelem++ )
00146 {
00147 if( abund.lgAbunTabl[nelem] )
00148 {
00149 abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)*
00150 dense.gas_phase[ipHYDROGEN];
00151
00152 hold = abun/dense.gas_phase[nelem];
00153 dense.gas_phase[nelem] = (realnum)abun;
00154
00155 for( ion=0; ion < (nelem + 2); ion++ )
00156 {
00157 dense.xIonDense[nelem][ion] *= (realnum)hold;
00158 }
00159 }
00160 }
00161 }
00162
00163
00164
00165
00166 if( !dense.lgDenFlucOn )
00167 {
00168
00169 lgChange = true;
00170 if( nzone <= 1 )
00171 {
00172 OldAbun = 1.;
00173 FacAbun = 1.;
00174 if( dense.lgDenFlucRadius )
00175 {
00176
00177 FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
00178 dense.flcPhase) + dense.csecnd;
00179 }
00180 else
00181 {
00182
00183 FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
00184 dense.flcPhase) + dense.csecnd;
00185 }
00186 }
00187 else
00188 {
00189 OldAbun = FacAbunSav;
00190
00191 if( dense.lgDenFlucRadius )
00192 {
00193
00194 FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
00195 dense.flcPhase) + dense.csecnd;
00196 }
00197 else
00198 {
00199
00200 FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
00201 dense.flcPhase) + dense.csecnd;
00202 }
00203 FacAbun = FacAbunSav/OldAbun;
00204 }
00205 }
00206 else
00207 {
00208
00209 FacAbun = 1.;
00210 }
00211
00212
00213
00214 if( lgChange )
00215 {
00216
00217
00218 for( nelem=ipHYDROGEN; nelem <= ipHELIUM; ++nelem )
00219 {
00220 dense.xMolecules[nelem] *= (realnum)pressure_change_factor;
00221 dense.gas_phase[nelem] *= (realnum)pressure_change_factor;
00222 for( ion=0; ion < (nelem + 2); ion++ )
00223 {
00224
00225 dense.xIonDense[nelem][ion] *= (realnum)pressure_change_factor;
00226 }
00227 }
00228
00229
00230
00231 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
00232 {
00233 dense.xMolecules[nelem] *= (realnum)(pressure_change_factor*FacAbun);
00234 dense.gas_phase[nelem] *= (realnum)(pressure_change_factor*FacAbun);
00235 for( ion=0; ion < (nelem + 2); ion++ )
00236 {
00237
00238 dense.xIonDense[nelem][ion] *= (realnum)(pressure_change_factor*FacAbun);
00239 }
00240 }
00241
00242 dense.eden *= pressure_change_factor;
00243
00244
00245
00246 TempChange(phycon.te , false);
00247
00248
00249 for(mol = 0; mol < N_H_MOLEC; mol++)
00250 {
00251 hmi.Hmolec[mol] *= (realnum)pressure_change_factor;
00252 }
00253 hmi.H2_total *= (realnum)pressure_change_factor;
00254
00255
00256
00257
00258 for( ion=0; ion < mole.num_comole_calc; ion++ )
00259 {
00260 COmole[ion]->hevmol *= (realnum)(pressure_change_factor*FacAbun);
00261
00262
00263 ASSERT( !isnan( COmole[ion]->hevmol ) );
00264 }
00265 }
00266
00267 if( trace.lgTrace )
00268 {
00269 fprintf( ioQQQ,
00270 " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
00271 hold, dense.gas_phase[ipHYDROGEN], geometry.FillFac );
00272
00273 if( trace.lgNeBug )
00274 {
00275 fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n",
00276 edensave, dense.eden, edensave/dense.eden );
00277 }
00278 }
00279
00280 {
00281
00282 enum {DEBUG_LOC=false};
00283
00284 if( DEBUG_LOC && (nzone>215) )
00285 {
00286 fprintf( ioQQQ,
00287 "%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n",
00288 radius.depth,
00289 pressure.PresTotlCurr,
00290 pressure.PresTotlInit + pressure.PresInteg,
00291 pressure.PresTotlInit,
00292 pressure.PresGasCurr,
00293 pressure.PresRamCurr,
00294 pressure.pres_radiation_lines_curr,
00295
00296 pressure.PresInteg - pressure.pinzon,
00297 wind.windv/1e5,
00298 sqrt(5.*pressure.PresGasCurr/3./dense.xMassDensity)/1e5,
00299 TorF(conv.lgConvPres) );
00300 }
00301 }
00302
00303 return lgChange;
00304 }
00305
00306
00307
00308
00309 bool lgConvPres(void)
00310 {
00311 double dnew,
00312 term;
00313 bool lgRet;
00314
00315 DEBUG_ENTRY( "lgConvPres()" );
00316
00317
00318
00319 pressure.PresTotlCorrect = 0.;
00320
00321
00322
00323
00324 if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00325 {
00326
00327 if( radius.glbdst < 0. )
00328 {
00329 fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
00330 fprintf( ioQQQ, " This is routine lgConvPres, GLBDST is%10.2e\n",
00331 radius.glbdst );
00332 cdEXIT(EXIT_FAILURE);
00333 }
00334 pressure_change_factor = (radius.glbden*pow(radius.glbrad/(radius.glbdst),radius.glbpow))/
00335 dense.gas_phase[ipHYDROGEN];
00336 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00337 }
00338 else if( cosmology.lgDo )
00339 {
00340
00341 dnew = GetDensity( cosmology.redshift_current );
00342 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
00343 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00344 }
00345 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) ) {
00346
00347
00348 if ( wind.lgStatic() )
00349 {
00350
00351 fprintf( ioQQQ, " PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" );
00352
00353 TotalInsanity();
00354 }
00355
00356
00357 else if( wind.lgBallistic() )
00358 {
00359
00360
00361
00362
00363 if( nzone > 1 )
00364 {
00365
00366
00367 if( trace.lgTrace && trace.lgWind )
00368 {
00369 fprintf(ioQQQ," lgConvPres sets AccelGravity %.3e lgDisk?%c\n",
00370 wind.AccelGravity ,
00371 TorF(wind.lgDisk) );
00372 }
00373
00374
00375
00376
00377
00378 term = POW2(struc.windv[nzone-2]) + 2.*(wind.AccelTotalOutward - wind.AccelGravity)* radius.drad;
00379
00380
00381 if( term <= 1e3 )
00382 {
00383
00384
00385 wind.lgVelPos = false;
00386 }
00387 else
00388 {
00389
00390
00391
00392
00393
00394
00395
00396 term = sqrt(term);
00397 if (wind.windv > 0)
00398 {
00399 wind.windv = (realnum) term;
00400 }
00401 else
00402 {
00403 wind.windv = -(realnum) term;
00404 }
00405 wind.lgVelPos = true;
00406 }
00407
00408 if( trace.lgTrace && trace.lgWind )
00409 {
00410 fprintf(ioQQQ," lgConvPres new wind V zn%li %.3e AccelTotalOutward %.3e AccelGravity %.3e\n",
00411 nzone,wind.windv, wind.AccelTotalOutward, wind.AccelGravity );
00412 }
00413 }
00414 else
00415 {
00416 pressure_change_factor = 1.;
00417 }
00418
00419
00420 pressure_change_factor = wind.emdot/(wind.windv*dense.gas_phase[ipHYDROGEN])/radius.r1r0sq;
00421 pressure.PresTotlCorrect = pressure.PresTotlCurr * pressure_change_factor;
00422 }
00423
00424
00425 else
00426 {
00427
00428 pressure_change_factor = DynaPresChngFactor();
00429 }
00430 }
00431
00432 else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
00433 {
00434
00435 if( dense.lgDenFlucRadius )
00436 {
00437 pressure_change_factor = (dense.cfirst*cos(radius.depth*dense.flong+dense.flcPhase) +
00438 dense.csecnd)/dense.gas_phase[ipHYDROGEN];
00439 }
00440 else
00441 {
00442 pressure_change_factor = (dense.cfirst*cos(colden.colden[ipCOL_HTOT]*dense.flong+dense.flcPhase) +
00443 dense.csecnd)/dense.gas_phase[ipHYDROGEN];
00444 }
00445 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00446 }
00447
00448 else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
00449 {
00450
00451 dnew = dense.den0*pow(radius.Radius/radius.rinner,(double)dense.DensityPower);
00452 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
00453 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00454 }
00455
00456 else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
00457 {
00458
00459 dnew = dense.den0*pow(1. + radius.depth/dense.rscale,(double)dense.DensityPower);
00460 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
00461 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00462 }
00463
00464 else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
00465 {
00466
00467 dnew = dense.den0*pow(1.f + colden.colden[ipCOL_HTOT]/
00468 dense.rscale,dense.DensityPower);
00469 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
00470 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00471 }
00472
00473 else if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
00474 {
00475
00476 if( pressure.lgContRadPresOn )
00477 {
00478
00479
00480 pressure.PresTotlCorrect = pressure.PresTotlInit + pressure.PresInteg;
00481 }
00482 else
00483 {
00484
00485 pressure.PresTotlCorrect = pressure.PresTotlInit*
00486
00487 pow(radius.Radius/radius.rinner,(double)pressure.PresPowerlaw);
00488 }
00489
00490 if( dark.lgNFW_Set || pressure.gravity_symmetry >=0 )
00491 {
00492 fixit();
00493
00494 pressure.PresTotlCorrect += pressure.IntegRhoGravity;
00495 }
00496
00497
00498 pressure_change_factor = pressure.PresTotlCorrect / pressure.PresTotlCurr;
00499 }
00500
00501 else if( strncmp( dense.chDenseLaw ,"DLW" , 3) == 0 )
00502 {
00503 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00504 {
00505
00506 pressure_change_factor = dense_fabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN];
00507 }
00508 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00509 {
00510
00511
00512 pressure_change_factor = dense_tabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN];
00513 }
00514 else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
00515 {
00516
00517 pressure_change_factor = dense_parametric_wind(radius.Radius)/dense.gas_phase[ipHYDROGEN];
00518 }
00519 else
00520 {
00521 fprintf( ioQQQ, " Insanity, lgConvPres gets chCPres=%4.4s\n",
00522 dense.chDenseLaw );
00523 cdEXIT(EXIT_FAILURE);
00524 }
00525 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00526
00527
00528
00529 if( pressure_change_factor > 3. || pressure_change_factor< 1./3 )
00530 {
00531 static bool lgWARN2BIG=false;
00532 if( !lgWARN2BIG )
00533 {
00534 lgWARN2BIG = true;
00535 fprintf(ioQQQ,"\n\n >========== Warning! The tabulated or functional change in density as a function of depth was VERY large. This is zone %li.\n",nzone);
00536 fprintf(ioQQQ," >========== Warning! This will cause convergence problems. \n");
00537 fprintf(ioQQQ," >========== Warning! The current radius is %.3e. \n",radius.Radius);
00538 fprintf(ioQQQ," >========== Warning! The current depth is %.3e. \n",radius.depth);
00539 fprintf(ioQQQ," >========== Warning! Consider using a more modest change in density vs radius. \n\n\n");
00540 }
00541 }
00542 }
00543
00544 else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
00545 {
00546
00547 pressure_change_factor = 1.;
00548 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00549 }
00550
00551 else
00552 {
00553 fprintf( ioQQQ, " Unknown pressure law=%s= This is a major internal error.\n",
00554 dense.chDenseLaw );
00555 ShowMe();
00556 cdEXIT(EXIT_FAILURE);
00557 }
00558
00559
00560
00561 ASSERT( pressure.PresTotlCorrect > FLT_MIN );
00562
00563
00564 if( pressure_change_factor > 1. + conv.PressureErrorAllowed || pressure_change_factor < 1. - conv.PressureErrorAllowed )
00565 {
00566 lgRet = false;
00567 conv.lgConvPres = false;
00568 }
00569 else
00570 {
00571 lgRet = true;
00572 conv.lgConvPres = true;
00573 }
00574
00575 return lgRet;
00576 }