00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "cddrive.h"
00007 #include "physconst.h"
00008 #include "iso.h"
00009 #include "taulines.h"
00010 #include "hydrogenic.h"
00011 #include "struc.h"
00012 #include "dynamics.h"
00013 #include "prt.h"
00014 #include "hyperfine.h"
00015 #include "dense.h"
00016 #include "magnetic.h"
00017 #include "continuum.h"
00018 #include "geometry.h"
00019 #include "h2.h"
00020 #include "co.h"
00021 #include "he.h"
00022 #include "grains.h"
00023 #include "atomfeii.h"
00024 #include "pressure.h"
00025 #include "stopcalc.h"
00026 #include "conv.h"
00027 #include "mean.h"
00028 #include "ca.h"
00029 #include "thermal.h"
00030 #include "atoms.h"
00031 #include "wind.h"
00032 #include "opacity.h"
00033 #include "timesc.h"
00034 #include "trace.h"
00035 #include "colden.h"
00036 #include "secondaries.h"
00037 #include "hmi.h"
00038 #include "radius.h"
00039 #include "phycon.h"
00040 #include "called.h"
00041 #include "mole.h"
00042 #include "rfield.h"
00043 #include "ionbal.h"
00044 #include "atmdat.h"
00045 #include "lines.h"
00046 #include "molcol.h"
00047 #include "input.h"
00048 #include "rt.h"
00049 #include "iterations.h"
00050 #include "cosmology.h"
00051 #include "deuterium.h"
00052
00053
00054 static double h2plus_heat_save, HeatH2Dish_used_save, HeatH2Dexc_used_save,
00055 hmihet_save, hmitot_save, H2_Solomon_dissoc_rate_used_H2g_save,deriv_HeatH2Dexc_used_save,
00056 H2_Solomon_dissoc_rate_used_H2s_save, H2_H2g_to_H2s_rate_used_save, H2_photodissoc_used_H2s_save,
00057 H2_photodissoc_used_H2g_save,
00058 UV_Cont_rel2_Draine_DB96_face,
00059 UV_Cont_rel2_Draine_DB96_depth , UV_Cont_rel2_Habing_TH85_face ,
00060 **saveMoleSource , **saveMoleSink;
00061 static realnum ***SaveMoleChTrRate;
00062
00063
00064 static realnum xIonFsave[LIMELM][LIMELM+1];
00065 static realnum deutDenseSave0, deutDenseSave1;
00066 static double HeatSave[LIMELM][LIMELM];
00067 static realnum supsav[LIMELM][LIMELM],p2nit,d5200r;
00068
00069 static double drSave , drNextSave;
00070
00071
00072 static long int IonLowSave[LIMELM],
00073 IonHighSave[LIMELM];
00074
00075
00076
00077 static bool lgHNSAV = false;
00078
00079 static double ***HOpacRatSav;
00080
00081 static double ortho_save , para_save;
00082 static double ***hnsav,
00083 edsav;
00084
00085 static realnum gas_phase_save[LIMELM];
00086 static double *den_save;
00087
00088
00089 void IterStart(void)
00090 {
00091 long int i,
00092 ion,
00093 ion2,
00094 ipISO,
00095 n ,
00096 nelem;
00097 double fhe1nx,
00098 ratio;
00099
00100 DEBUG_ENTRY( "IterStart()" );
00101
00102
00103
00104
00105
00106 if( !lgHNSAV )
00107 {
00108
00109 lgHNSAV = true;
00110
00111 HOpacRatSav = (double ***)MALLOC(sizeof(double **)*NISO );
00112
00113 hnsav = (double ***)MALLOC(sizeof(double **)*NISO );
00114
00115 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00116 {
00117
00118 HOpacRatSav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM );
00119
00120 hnsav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM );
00121
00122
00123 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00124 {
00125 HOpacRatSav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso_sp[ipISO][nelem].numLevels_max+1) );
00126
00127 hnsav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso_sp[ipISO][nelem].numLevels_max+1) );
00128 }
00129 }
00130 saveMoleSource =
00131 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
00132 saveMoleSink =
00133 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
00134 SaveMoleChTrRate =
00135 (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM );
00136 for(nelem=0; nelem<LIMELM; ++nelem )
00137 {
00138
00139 saveMoleSource[nelem] =
00140 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
00141 saveMoleSink[nelem] =
00142 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
00143 SaveMoleChTrRate[nelem] =
00144 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) );
00145 for( ion=0; ion<nelem+2; ++ion )
00146 {
00147 SaveMoleChTrRate[nelem][ion] =
00148 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) );
00149 }
00150 }
00151 den_save = (double*)MALLOC(sizeof(double)*(mole_global.num_calc) );
00152 }
00153
00154
00155 if( trace.lgTrace )
00156 {
00157 fprintf( ioQQQ, " IterStart called.\n" );
00158 }
00159
00160
00161 if( iteration > iterations.itermx )
00162 iterations.lgLastIt = true;
00163
00164
00165 rt.lgMaserCapHit = false;
00166
00167 rt.lgMaserSetDR = false;
00168
00169
00170 atmdat.HCharHeatMax = 0.;
00171 atmdat.HCharCoolMax = 0.;
00172
00173
00174 for(nelem=0; nelem<HS_NZ; ++nelem )
00175 {
00176 atmdat.lgHCaseBOK[0][nelem] = true;
00177 atmdat.lgHCaseBOK[1][nelem] = true;
00178 }
00179
00180
00181 for( long ipSpecies=0; ipSpecies<nSpecies; ++ipSpecies )
00182 {
00183 dBaseSpecies[ipSpecies].numLevels_local = dBaseSpecies[ipSpecies].numLevels_max;
00184 dBaseSpecies[ipSpecies].CoolTotal = 0.;
00185 }
00186
00187
00188 strncpy( StopCalc.chReasonStop, "reason not specified.",
00189 sizeof(StopCalc.chReasonStop) );
00190
00191
00192 he.nzone = 0;
00193 he.frac_he0dest_23S = 0.;
00194 he.frac_he0dest_23S_photo = 0.;
00195
00196 dense.EdenMax = 0.;
00197 dense.EdenMin = BIGFLOAT;
00198
00199 timesc.time_H2_Dest_longest = 0.;
00200 timesc.time_H2_Form_longest = 0.;
00201
00202 timesc.time_H2_Dest_here = -1.;
00203 timesc.time_H2_Form_here = 0.;
00204
00205 timesc.BigCOMoleForm = 0.;
00206 timesc.TimeH21cm = 0.;
00207 thermal.CoolHeatMax = 0.;
00208 hydro.HCollIonMax = 0.;
00209
00210 atmdat.HIonFracMax = 0.;
00211
00212
00213 hydro.pop2mx = 0.;
00214 hydro.lgHiPop2 = false;
00215
00216 hydro.nLyaHot = 0;
00217 hydro.TLyaMax = 0.;
00218
00219
00220
00221 pressure.PresInteg = 0.;
00222 pressure.pinzon = 0.;
00223 pressure.PresIntegElecThin = 0.;
00224 PresTotCurrent();
00225
00226 dynamics.HeatMax = 0.;
00227 dynamics.CoolMax = 0.;
00228 if( dynamics.lgAdvection )
00229 DynaIterStart();
00230
00231
00232 pressure.pbeta = 0.;
00233 pressure.RadBetaMax = 0.;
00234 pressure.lgPradCap = false;
00235 pressure.lgPradDen = false;
00236
00237 pressure.lgSonicPoint = false;
00238 pressure.lgStrongDLimbo = false;
00239
00240
00241 timesc.tcmptn = 1.e0/(rfield.qtot*6.65e-25*dense.eden);
00242 timesc.time_therm_long = 1.5*dense.pden*BOLTZMANN*phycon.te/thermal.ctot;
00243
00244
00245 dense.xMassTotal = 0.;
00246
00247 for( nelem=0; nelem < LIMELM; nelem++ )
00248 {
00249
00250 if( dense.lgElmtOn[nelem] )
00251 {
00252
00253
00254 for( ion=0; ion < (nelem + 2); ion++ )
00255 {
00256 xIonFsave[nelem][ion] = dense.xIonDense[nelem][ion];
00257 }
00258
00259 for( ion=0; ion < LIMELM; ion++ )
00260 {
00261 HeatSave[nelem][ion] = thermal.heating[nelem][ion];
00262 }
00263 }
00264 }
00265
00266 deutDenseSave0 = deut.xIonDense[0];
00267 deutDenseSave1 = deut.xIonDense[1];
00268
00269
00270 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00271 {
00272
00273 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00274 {
00275 if( dense.lgElmtOn[nelem] )
00276 {
00277 for( n=0; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
00278 {
00279 HOpacRatSav[ipISO][nelem][n] = iso_sp[ipISO][nelem].fb[n].ConOpacRatio;
00280 hnsav[ipISO][nelem][n] = iso_sp[ipISO][nelem].st[n].Pop();
00281 }
00282 }
00283 for( vector<two_photon>::iterator tnu = iso_sp[ipISO][nelem].TwoNu.begin(); tnu != iso_sp[ipISO][nelem].TwoNu.end(); ++tnu )
00284 tnu->induc_dn_max = 0.;
00285 }
00286 }
00287
00288 rfield.ipEnergyBremsThin = 0;
00289 rfield.EnergyBremsThin = 0.;
00290
00291 p2nit = atoms.p2nit;
00292 d5200r = atoms.d5200r;
00293 atoms.nNegOI = 0;
00294 for( i=0; i< N_OI_LEVELS; ++i )
00295 atoms.popoi[i] = 0.;
00296
00297
00298 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00299 {
00300
00301 gas_phase_save[nelem] = dense.gas_phase[nelem];
00302 IonLowSave[nelem] = dense.IonLow[nelem];
00303 IonHighSave[nelem] = dense.IonHigh[nelem];
00304 }
00305
00306
00307
00308 edsav = dense.eden;
00309
00310
00311 for( i=0; i<mole_global.num_calc; i++)
00312 {
00313 den_save[i] = mole.species[i].den;
00314 }
00315 ortho_save = h2.ortho_density;
00316 para_save = h2.para_density;
00317
00318 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
00319 {
00320
00321 for( ion=0; ion<nelem+2; ++ion )
00322 {
00323
00324 saveMoleSource[nelem][ion] = mole.source[nelem][ion];
00325 saveMoleSink[nelem][ion] = mole.sink[nelem][ion];
00326
00327
00328
00329 for( ion2=0; ion2<nelem+2; ++ion2 )
00330 {
00331 SaveMoleChTrRate[nelem][ion][ion2] = mole.xMoleChTrRate[nelem][ion][ion2];
00332 }
00333 }
00334 }
00335
00336 hmihet_save = hmi.hmihet;
00337 hmitot_save = hmi.hmitot;
00338
00339 h2plus_heat_save = hmi.h2plus_heat;
00340
00341 HeatH2Dish_used_save = hmi.HeatH2Dish_used;
00342 HeatH2Dexc_used_save = hmi.HeatH2Dexc_used;
00343
00344 deriv_HeatH2Dexc_used_save = hmi.deriv_HeatH2Dexc_used;
00345 H2_Solomon_dissoc_rate_used_H2g_save = hmi.H2_Solomon_dissoc_rate_used_H2g;
00346 H2_Solomon_dissoc_rate_used_H2s_save = hmi.H2_Solomon_dissoc_rate_used_H2s;
00347
00348 H2_H2g_to_H2s_rate_used_save = hmi.H2_H2g_to_H2s_rate_used;
00349
00350 H2_photodissoc_used_H2s_save = hmi.H2_photodissoc_used_H2s;
00351 H2_photodissoc_used_H2g_save = hmi.H2_photodissoc_used_H2g;
00352
00353 UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face;
00354 UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_depth;
00355 UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_face;
00356
00357
00358 drSave = radius.drad;
00359 drNextSave = radius.drNext;
00360
00361 radius.dr_min_last_iter = BIGFLOAT;
00362 radius.dr_max_last_iter = 0.;
00363
00364 geometry.nprint = iterations.IterPrnt[iteration-1];
00365
00366 colden.TotMassColl = 0.;
00367 colden.tmas = 0.;
00368 colden.wmas = 0.;
00369 colden.rjnmin = FLT_MAX;
00370 colden.ajmmin = FLT_MAX;
00371
00372 thermal.nUnstable = 0;
00373 thermal.lgUnstable = false;
00374
00375 rfield.extin_mag_B_point = 0.;
00376 rfield.extin_mag_V_point = 0.;
00377 rfield.extin_mag_B_extended = 0.;
00378 rfield.extin_mag_V_extended = 0.;
00379
00380
00381 for( i=0; i < rfield.nupper+1; i++ )
00382 {
00383
00384
00385 rfield.SummedDifSave[i] = rfield.SummedDif[i];
00386 rfield.ConEmitReflec[0][i] = 0.;
00387 rfield.ConEmitOut[0][i] = 0.;
00388 rfield.ConInterOut[i] = 0.;
00389
00390 rfield.otssav[i][0] = rfield.otscon[i];
00391 rfield.otssav[i][1] = rfield.otslin[i];
00392 rfield.outlin[0][i] = 0.;
00393 rfield.outlin_noplot[i] = 0.;
00394 rfield.ConOTS_local_OTS_rate[i] = 0.;
00395 rfield.ConOTS_local_photons[i] = 0.;
00396 rfield.DiffuseEscape[i] = 0.;
00397
00398
00399 opac.opacity_abs_savzon1[i] = opac.opacity_abs[i];
00400 opac.opacity_sct_savzon1[i] = opac.opacity_sct[i];
00401
00402
00403 opac.TauAbsFace[i] = opac.taumin;
00404 opac.TauScatFace[i] = opac.taumin;
00405
00406
00407
00408
00409
00410
00411
00412
00413 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin);
00414
00415
00416 opac.ExpmTau[i] = (realnum)opac.ExpZone[i];
00417
00418
00419
00420 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i] );
00421 }
00422
00423 t_fe2ovr_la::Inst().zero_opacity();
00424
00425
00426
00427
00428 mean.MeanZero();
00429
00430
00431 for( i=0; i < NCOLD; i++ )
00432 {
00433 colden.colden[i] = 0.;
00434 }
00435 colden.He123S = 0.;
00436 colden.H0_ov_Tspin = 0.;
00437 colden.OH_ov_Tspin = 0.;
00438 colden.coldenH2_ov_vel = 0.;
00439
00440
00441 colden.H0_21cm_upper =0;
00442 colden.H0_21cm_lower =0;
00443
00444 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00445 {
00446 (*diatom)->ortho_colden = 0.;
00447 (*diatom)->para_colden = 0.;
00448 }
00449
00450 for( i=0; i < mole_global.num_calc; i++ )
00451 {
00452 mole.species[i].column = 0.;
00453 }
00454 for( i=0; i < 5; i++ )
00455 {
00456 colden.C2Pops[i] = 0.;
00457 colden.C2Colden[i] = 0.;
00458
00459 colden.Si2Pops[i] = 0.;
00460 colden.Si2Colden[i] = 0.;
00461 }
00462 for( i=0; i < 3; i++ )
00463 {
00464
00465 colden.C1Pops[i] = 0.;
00466 colden.C1Colden[i] = 0.;
00467
00468 colden.O1Pops[i] = 0.;
00469 colden.O1Colden[i] = 0.;
00470
00471 colden.C3Pops[i] = 0.;
00472 }
00473 for( i=0; i < 4; i++ )
00474 {
00475
00476 colden.C3Colden[i] = 0.;
00477 }
00478
00479
00480 colden.dlnenp = 0.;
00481 colden.dlnenHep = 0.;
00482 colden.dlnenHepp = 0.;
00483 colden.dlnenCp = 0.;
00484 colden.H0_ov_Tspin = 0.;
00485 colden.OH_ov_Tspin = 0.;
00486
00487
00488 for( unsigned i = 0; i < mole.species.size(); ++i )
00489 {
00490 if( mole.species[i].levels != NULL )
00491 {
00492 for( qList::iterator st = mole.species[i].levels->begin(); st != mole.species[i].levels->end(); ++st )
00493 {
00494 (*st).ColDen() = 0.;
00495 }
00496 }
00497 }
00498
00499
00500 molcol("ZERO",ioQQQ);
00501
00502
00503 thermal.FreeFreeTotHeat = 0.;
00504
00505 thermal.thist = 0.;
00506 thermal.tlowst = 1e20f;
00507
00508 wind.AccelAver = 0.;
00509 wind.acldr = 0.;
00510 ionbal.ilt = 0;
00511 ionbal.iltln = 0;
00512 ionbal.ilthn = 0;
00513 ionbal.ihthn = 0;
00514 ionbal.ifail = 0;
00515
00516 secondaries.SecHIonMax = 0.;
00517 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00518 {
00519 for( ion=0; ion<nelem+1; ++ion )
00520 {
00521 supsav[nelem][ion] = secondaries.csupra[nelem][ion];
00522 }
00523 }
00524 secondaries.x12sav = secondaries.x12tot;
00525 secondaries.hetsav = secondaries.HeatEfficPrimary;
00526 secondaries.savefi = secondaries.SecIon2PrimaryErg;
00527 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00528 {
00529 for( ion=0; ion<nelem+1; ++ion )
00530 {
00531 ionbal.CompRecoilHeatRateSave[nelem][ion] = ionbal.CompRecoilHeatRate[nelem][ion];
00532 ionbal.CompRecoilIonRateSave[nelem][ion] = ionbal.CompRecoilIonRate[nelem][ion];
00533 }
00534 }
00535
00536
00537 conv.nTeFail = 0;
00538 conv.nTotalFailures = 0;
00539 conv.nPreFail = 0;
00540 conv.failmx = 0.;
00541 conv.nIonFail = 0;
00542 conv.nPopFail = 0;
00543 conv.nNeFail = 0;
00544 conv.nGrainFail = 0;
00545 conv.nChemFail = 0;
00546 conv.dCmHdT = 0.;
00547
00548
00549 lgAbort = false;
00550
00551 GrainStartIter();
00552
00553 rfield.comtot = 0.;
00554
00555 co.codfrc = 0.;
00556 co.codtot = 0.;
00557
00558 hmi.HeatH2DexcMax = 0.;
00559 hmi.CoolH2DexcMax = 0.;
00560 hmi.h2dfrc = 0.;
00561 hmi.h2line_cool_frac = 0.;
00562 hmi.h2dtot = 0.;
00563 thermal.HeatLineMax = 0.;
00564 thermal.GBarMax = 0.;
00565 hyperfine.cooling_max = 0.;
00566 hydro.cintot = 0.;
00567 geometry.lgZoneTrp = false;
00568
00569 hmi.h2pmax = 0.;
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580 ASSERT( LineSave.nsum > 0);
00581 ASSERT( LineSave.nsumAllocated == LineSave.nsum );
00582
00583
00584 for( i=0; i < LineSave.nsum; i++ )
00585 {
00586 LineSv[i].SumLine[0] = 0.;
00587 LineSv[i].SumLine[1] = 0.;
00588 LineSv[i].emslin[0] = 0.;
00589 LineSv[i].emslin[1] = 0.;
00590 }
00591
00592
00593 hydro.cintot = 0.;
00594 thermal.totcol = 0.;
00595 rfield.comtot = 0.;
00596 thermal.FreeFreeTotHeat = 0.;
00597 thermal.power = 0.;
00598 thermal.HeatLineMax = 0.;
00599 thermal.GBarMax = 0.;
00600 hyperfine.cooling_max = 0.;
00601
00602 hmi.h2pmax = 0.;
00603
00604 co.codfrc = 0.;
00605 co.codtot = 0.;
00606
00607 hmi.h2dfrc = 0.;
00608 hmi.h2line_cool_frac = 0.;
00609 hmi.h2dtot = 0.;
00610 timesc.sound = 0.;
00611
00612 if( LineSave.lgNormSet )
00613 {
00614
00615
00616 LineSave.ipNormWavL = 0;
00617 LineSave.ipNormWavL = cdLine( LineSave.chNormLab , LineSave.WavLNorm , &fhe1nx , &ratio );
00618 if( LineSave.ipNormWavL < 0 )
00619 {
00620
00621 fprintf( ioQQQ, "PROBLEM could not find the normalisation line.\n");
00622 fprintf( ioQQQ,
00623 "IterStart could not find a line with a wavelength of %f and a label of %s\n",
00624 LineSave.WavLNorm, LineSave.chNormLab );
00625 fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
00626 fprintf( ioQQQ, "Sorry.\n");
00627 cdEXIT(EXIT_FAILURE);
00628 }
00629 }
00630 else
00631 {
00632
00633
00634 LineSave.ipNormWavL = 0;
00635 LineSave.ipNormWavL = cdLine( "TOTL" , 4861.36f , &fhe1nx , &ratio );
00636 if( LineSave.ipNormWavL < 0 )
00637 {
00638
00639 fprintf( ioQQQ, "PROBLEM could not find the default normalisation line.\n");
00640 fprintf( ioQQQ,
00641 "IterStart could not find a line with a wavelength of 4861 and a label of TOTL\n" );
00642 fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
00643 fprintf( ioQQQ, "Sorry.\n");
00644 cdEXIT(EXIT_FAILURE);
00645 }
00646 }
00647
00648
00649
00650
00651 if( iteration == 1 && StopCalc.nstpl )
00652 {
00653
00654 for( long int nStopLine=0; nStopLine < StopCalc.nstpl; nStopLine++ )
00655 {
00656 double relint, absint ;
00657
00658
00659 StopCalc.ipStopLin1[nStopLine] = cdLine( StopCalc.chStopLabel1[nStopLine],
00660
00661 StopCalc.StopLineWl1[nStopLine], &relint, &absint );
00662
00663 if( StopCalc.ipStopLin1[nStopLine]<0 )
00664 {
00665 fprintf( ioQQQ,
00666 " IterStart could not find first line in STOP LINE command, number %ld with label *%s* and wl %.1f\n",
00667 StopCalc.ipStopLin1[nStopLine] ,
00668 StopCalc.chStopLabel1[nStopLine],
00669 StopCalc.StopLineWl1[nStopLine]);
00670 cdEXIT(EXIT_FAILURE);
00671 }
00672
00673 StopCalc.ipStopLin2[nStopLine] = cdLine( StopCalc.chStopLabel2[nStopLine],
00674
00675 StopCalc.StopLineWl2[nStopLine], &relint, &absint );
00676
00677 if( StopCalc.ipStopLin2[nStopLine] < 0 )
00678 {
00679 fprintf( ioQQQ,
00680 " IterStart could not find second line in STOP LINE command, line number %ld with label *%s* and wl %.1f\n",
00681 StopCalc.ipStopLin1[nStopLine] ,
00682 StopCalc.chStopLabel1[nStopLine],
00683 StopCalc.StopLineWl1[nStopLine]);
00684 cdEXIT(EXIT_FAILURE);
00685 }
00686
00687 if( trace.lgTrace )
00688 {
00689 fprintf( ioQQQ,
00690 " stop line 1 is number %5ld wavelength is %f label is %4.4s\n",
00691 StopCalc.ipStopLin1[nStopLine],
00692 LineSv[StopCalc.ipStopLin1[nStopLine]].wavelength,
00693 LineSv[StopCalc.ipStopLin1[nStopLine]].chALab );
00694 fprintf( ioQQQ,
00695 " stop line 2 is number %5ld wavelength is %f label is %4.4s\n",
00696 StopCalc.ipStopLin2[nStopLine],
00697 LineSv[StopCalc.ipStopLin2[nStopLine]].wavelength,
00698 LineSv[StopCalc.ipStopLin2[nStopLine]].chALab );
00699 }
00700 }
00701 }
00702
00703
00704 if( prt.lgPrtLastIt )
00705 {
00706 if( iteration == 1 )
00707 {
00708 called.lgTalk = false;
00709 }
00710
00711
00712
00713
00714
00715 if( iterations.lgLastIt && !prt.lgPrtStart && !called.lgTalkForcedOff )
00716 {
00717 called.lgTalk = called.lgTalkSave;
00718 }
00719 }
00720
00721 if( opac.lgCaseB )
00722 {
00723 if( trace.lgTrace )
00724 {
00725 fprintf( ioQQQ, " IterStart does not change mid-zone optical depth since CASE B\n" );
00726 }
00727 }
00728
00729
00730 hydro.FracInd = 0.;
00731 hydro.fbul = 0.;
00732
00733
00734
00735 for( i=ipH2p; i < (iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max - 1); i++ )
00736 {
00737 if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().Aul() <= iso_ctrl.SmallA )
00738 continue;
00739
00740 ratio = iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RecomInducRate*iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].PopLTE/
00741 SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RecomInducRate*iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].PopLTE +
00742 iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RadRecomb[ipRecRad]*
00743 iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RadRecomb[ipRecNetEsc]);
00744 if( ratio > hydro.FracInd )
00745 {
00746 hydro.FracInd = (realnum)ratio;
00747 hydro.ndclev = i;
00748 }
00749
00750 ratio = iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().pump()/
00751 (iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().pump() +
00752 iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().Aul());
00753
00754 if( ratio > hydro.fbul )
00755 {
00756 hydro.fbul = (realnum)ratio;
00757 hydro.nbul = i;
00758 }
00759 }
00760
00761 if( trace.lgTrace )
00762 fprintf( ioQQQ, " IterStart returns.\n" );
00763 return;
00764 }
00765
00766
00767
00768
00769
00770
00771 void IterRestart(void)
00772 {
00773 long int i,
00774 ion,
00775 ion2,
00776 ipISO ,
00777 n,
00778 nelem;
00779 double SumOTS;
00780
00781 DEBUG_ENTRY( "IterRestart()" );
00782
00783
00784
00785
00786
00787
00788
00789 if( StopCalc.TeFloor > 0. && !thermal.lgTemperatureConstantCommandParsed )
00790 {
00791 thermal.lgTemperatureConstant = false;
00792 thermal.ConstTemp = 0.;
00793 }
00794
00795 lgAbort = false;
00796
00797
00798 Magnetic_reinit();
00799
00800 opac.stimax[0] = 0.;
00801 opac.stimax[1] = 0.;
00802
00803 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00804 (*diatom)->H2_Reset();
00805
00806 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
00807 {
00808
00809 for( ion=0; ion<nelem+2; ++ion )
00810 {
00811
00812 mole.source[nelem][ion] = saveMoleSource[nelem][ion];
00813 mole.sink[nelem][ion] = saveMoleSink[nelem][ion];
00814
00815
00816
00817 for( ion2=0; ion2<nelem+2; ++ion2 )
00818 {
00819 mole.xMoleChTrRate[nelem][ion][ion2] = SaveMoleChTrRate[nelem][ion][ion2];
00820 }
00821 }
00822 }
00823
00824
00825 for( i=0; i<mole_global.num_calc; i++)
00826 {
00827 mole.species[i].den = den_save[i];
00828 }
00829 dense.updateXMolecules();
00830 hmi.H2_total = findspecieslocal("H2")->den + findspecieslocal("H2*")->den;
00831 hmi.HD_total = findspecieslocal("HD")->den;
00832
00833 h2.ortho_density = ortho_save;
00834 h2.para_density = para_save;
00835 {
00836 hmi.H2_total_f = (realnum)hmi.H2_total;
00837 h2.ortho_density_f = (realnum)h2.ortho_density;
00838 h2.para_density_f = (realnum)h2.para_density;
00839 }
00840
00841
00842 for( i=0; i < NCOLD; i++ )
00843 {
00844 colden.colden[i] = 0.;
00845 }
00846 colden.He123S = 0.;
00847 colden.coldenH2_ov_vel = 0.;
00848
00849 for( i=0; i < mole_global.num_calc; i++ )
00850 {
00851
00852 mole.species[i].column = 0.;
00853
00854 mole.species[i].xFracLim = 0.;
00855 mole.species[i].nAtomLim = -1;
00856 }
00857
00858
00859
00860 if( dynamics.lgAdvection )
00861 DynaIterEnd();
00862
00863 rfield.extin_mag_B_point = 0.;
00864 rfield.extin_mag_V_point = 0.;
00865 rfield.extin_mag_B_extended = 0.;
00866 rfield.extin_mag_V_extended = 0.;
00867
00868 hmi.hmihet = hmihet_save;
00869 hmi.hmitot = hmitot_save;
00870
00871 hmi.h2plus_heat = h2plus_heat_save;
00872 hmi.HeatH2Dish_used = HeatH2Dish_used_save;
00873 hmi.HeatH2Dexc_used = HeatH2Dexc_used_save;
00874
00875 hmi.deriv_HeatH2Dexc_used = (realnum)deriv_HeatH2Dexc_used_save;
00876 hmi.H2_Solomon_dissoc_rate_used_H2g = H2_Solomon_dissoc_rate_used_H2g_save;
00877 hmi.H2_Solomon_dissoc_rate_used_H2s = H2_Solomon_dissoc_rate_used_H2s_save;
00878 hmi.H2_H2g_to_H2s_rate_used = H2_H2g_to_H2s_rate_used_save;
00879 hmi.H2_photodissoc_used_H2s = H2_photodissoc_used_H2s_save;
00880 hmi.H2_photodissoc_used_H2g = H2_photodissoc_used_H2g_save;
00881
00882 hmi.UV_Cont_rel2_Draine_DB96_face = (realnum)UV_Cont_rel2_Draine_DB96_face;
00883 hmi.UV_Cont_rel2_Draine_DB96_depth = (realnum)UV_Cont_rel2_Draine_DB96_depth;
00884 hmi.UV_Cont_rel2_Habing_TH85_face = (realnum)UV_Cont_rel2_Habing_TH85_face;
00885
00886 rfield.ipEnergyBremsThin = 0;
00887 rfield.EnergyBremsThin = 0.;
00888 rfield.lgUSphON = false;
00889
00890 radius.glbdst = 0.;
00891
00892 radius.drad = drSave;
00893 radius.drad_mid_zone = drSave/2.;
00894 radius.drNext = drNextSave;
00895
00896
00897 radius.lgDrMinUsed = false;
00898
00899 continuum.cn4861 = continuum.sv4861;
00900 continuum.cn1216 = continuum.sv1216;
00901
00902
00903 continuum.totlsv = continuum.TotalLumin;
00904
00905
00906 if( (trace.nznbug == 0 && iteration >= trace.npsbug) && trace.lgTrOvrd )
00907 {
00908 if( trace.nTrConvg==0 )
00909 {
00910 trace.lgTrace = true;
00911 }
00912 else
00913
00914
00915 trace.nTrConvg = abs( trace.nTrConvg );
00916
00917 fprintf( ioQQQ, " IterRestart called.\n" );
00918 }
00919 else
00920 {
00921 trace.lgTrace = false;
00922 }
00923
00924
00925 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00926 {
00927 for( ion=0; ion<nelem+1; ++ion )
00928 {
00929 secondaries.csupra[nelem][ion] = supsav[nelem][ion];
00930 }
00931 }
00932 secondaries.x12tot = secondaries.x12sav;
00933 secondaries.HeatEfficPrimary = secondaries.hetsav;
00934 secondaries.SecIon2PrimaryErg = secondaries.savefi;
00935 for( nelem=0; nelem<LIMELM; ++nelem)
00936 {
00937 for( ion=0; ion<nelem+1; ++ion )
00938 {
00939 ionbal.CompRecoilHeatRate[nelem][ion] = ionbal.CompRecoilHeatRateSave[nelem][ion];
00940 ionbal.CompRecoilIonRate[nelem][ion] = ionbal.CompRecoilIonRateSave[nelem][ion];
00941 }
00942 }
00943
00944 wind.lgVelPos = true;
00945 wind.AccelMax = 0.;
00946 wind.AccelAver = 0.;
00947 wind.acldr = 0.;
00948 wind.windv = wind.windv0;
00949
00950 thermal.nUnstable = 0;
00951 thermal.lgUnstable = false;
00952
00953 pressure.pbeta = 0.;
00954 pressure.RadBetaMax = 0.;
00955 pressure.lgPradCap = false;
00956 pressure.lgPradDen = false;
00957
00958 pressure.lgSonicPoint = false;
00959 pressure.lgStrongDLimbo = false;
00960
00961 pressure.PresInteg = 0.;
00962 pressure.pinzon = 0.;
00963 pressure.PresIntegElecThin = 0.;
00964
00965 pressure.RhoGravity_dark = 0.;
00966 pressure.RhoGravity_self = 0.;
00967 pressure.RhoGravity_external = 0.;
00968 pressure.RhoGravity = 0.;
00969 pressure.IntegRhoGravity = 0.;
00970
00971 EdenChange( edsav );
00972 dense.EdenHCorr = edsav;
00973 dense.EdenHCorr_f = (realnum)dense.EdenHCorr;
00974 dense.EdenTrue = edsav;
00975
00976 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00977 {
00978
00979 dense.SetGasPhaseDensity( nelem, gas_phase_save[nelem] );
00980 dense.IonLow[nelem] = IonLowSave[nelem];
00981 dense.IonHigh[nelem] = IonHighSave[nelem];
00982 }
00983
00984
00985
00986 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00987 {
00988 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00989 {
00990 iso_sp[ipISO][nelem].Reset();
00991 }
00992 }
00993
00994 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00995 {
00996 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00997 {
00998 if( dense.lgElmtOn[nelem] )
00999 {
01000 for( n=ipH1s; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
01001 {
01002 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = HOpacRatSav[ipISO][nelem][n];
01003 iso_sp[ipISO][nelem].st[n].Pop() = hnsav[ipISO][nelem][n];
01004 }
01005 }
01006 }
01007 }
01008
01009 if( trace.lgTrace && trace.lgNeBug )
01010 {
01011 fprintf( ioQQQ, " EDEN set to%12.4e by IterRestart.\n",
01012 dense.eden );
01013 }
01014
01015 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
01016 {
01017 for( ion=0; ion < (nelem + 2); ion++ )
01018 {
01019 dense.xIonDense[nelem][ion] = xIonFsave[nelem][ion];
01020 }
01021 for( i=0; i < LIMELM; i++ )
01022 {
01023 thermal.heating[nelem][i] = HeatSave[nelem][i];
01024 }
01025 }
01026
01027 deut.xIonDense[0] = deutDenseSave0;
01028 deut.xIonDense[1] = deutDenseSave1;
01029
01030 GrainRestartIter();
01031
01032 rfield.resetCoarseTransCoef();
01033
01034
01035 for( i=0; i < rfield.nupper; i++ )
01036 {
01037
01038 rfield.flux_beam_const[i] = rfield.flux_beam_const_save[i];
01039
01040
01041 rfield.flux_beam_time[i] = rfield.flux_time_beam_save[i]*
01042 rfield.time_continuum_scale;
01043
01044 if( cosmology.lgDo )
01045 {
01046 double slope_ratio;
01047 double fac_old = TE1RYD*rfield.anu[i]/(CMB_TEMP*(1. + cosmology.redshift_start ));
01048 double fac_new = TE1RYD*rfield.anu[i]/(CMB_TEMP*(1. + cosmology.redshift_current));
01049
01050 if( fac_old > log10(DBL_MAX) )
01051 {
01052 slope_ratio = 0.;
01053 }
01054 else if( fac_new > log10(DBL_MAX) )
01055 {
01056 slope_ratio = 0.;
01057 }
01058 else if( fac_old > 1.e-5 )
01059 {
01060 slope_ratio = (exp(fac_old) - 1.)/(exp(fac_new) - 1.);
01061 }
01062 else
01063 {
01064 slope_ratio = (fac_old*(1. + fac_old/2.))/(fac_new*(1. + fac_new/2.));
01065 }
01066
01067 rfield.flux_isotropic[i] = rfield.flux_isotropic_save[i] *
01068 pow( ((realnum)1.f + cosmology.redshift_current)/((realnum)1.f + cosmology.redshift_start), (realnum)4.f ) * (realnum)slope_ratio;
01069 }
01070 else
01071 {
01072
01073 rfield.flux_isotropic[i] = rfield.flux_isotropic_save[i];
01074 }
01075
01076 rfield.flux[0][i] = rfield.flux_isotropic[i] + rfield.flux_beam_time[i] +
01077 rfield.flux_beam_const[i];
01078
01079 rfield.SummedDif[i] = rfield.SummedDifSave[i];
01080 rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
01081
01082 rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
01083 rfield.otscon[i] = rfield.otssav[i][0];
01084 rfield.otslin[i] = rfield.otssav[i][1];
01085 rfield.ConInterOut[i] = 0.;
01086 rfield.OccNumbBremsCont[i] = 0.;
01087 rfield.OccNumbDiffCont[i] = 0.;
01088 rfield.OccNumbContEmitOut[i] = 0.;
01089 rfield.outlin[0][i] = 0.;
01090 rfield.outlin_noplot[i] = 0.;
01091 rfield.ConOTS_local_OTS_rate[i] = 0.;
01092 rfield.ConOTS_local_photons[i] = 0.;
01093 rfield.DiffuseEscape[i] = 0.;
01094
01095 opac.opacity_abs[i] = opac.opacity_abs_savzon1[i];
01096 opac.OldOpacSave[i] = opac.opacity_abs_savzon1[i];
01097 opac.opacity_sct[i] = opac.opacity_sct_savzon1[i];
01098 opac.albedo[i] =
01099 opac.opacity_sct[i]/MAX2(1e-30,opac.opacity_sct[i] + opac.opacity_abs[i]);
01100 opac.tmn[i] = 1.;
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad/2.*geometry.DirectionalCosin);
01114
01115
01116 opac.ExpmTau[i] = (realnum)opac.ExpZone[i];
01117
01118
01119 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
01120 rfield.reflin[0][i] = 0.;
01121 rfield.ConEmitReflec[0][i] = 0.;
01122 rfield.ConEmitOut[0][i] = 0.;
01123 rfield.ConRefIncid[0][i] = 0.;
01124
01125
01126
01127 if( iteration > 1 )
01128 {
01129
01130 realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
01131 opac.E2TauAbsOut[i] = (realnum)e2( tau );
01132 ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
01133 }
01134 else
01135 opac.E2TauAbsOut[i] = 1.;
01136
01137 }
01138
01139
01140 RT_OTS_Update(&SumOTS);
01141
01142 thermal.FreeFreeTotHeat = 0.;
01143 atoms.p2nit = p2nit;
01144 atoms.d5200r = d5200r;
01145
01146 if( called.lgTalk )
01147 {
01148 fprintf( ioQQQ, "\f\n Start Iteration Number %ld %75.75s\n\n\n",
01149 iteration, input.chTitle );
01150 }
01151
01152
01153 FeIIReset();
01154 ASSERT(lgElemsConserved());
01155 return;
01156 }
01157
01158
01159 void IterEnd(void)
01160 {
01161
01162 DEBUG_ENTRY( "IterEnd()" );
01163
01164 if( lgAbort )
01165 {
01166 return;
01167 }
01168
01169
01170 double fac = radius.depth/radius.rinner;
01171 if( fac < 0.1 )
01172 {
01173 geometry.lgGeoPP = true;
01174 }
01175 else
01176 {
01177 geometry.lgGeoPP = false;
01178 }
01179
01180 if( iteration > dynamics.n_initial_relax && dynamics.lgTimeDependentStatic )
01181 {
01182
01183
01184
01185 double cumulative_factor = (dynamics.timestep/
01186 colden.TotMassColl);
01187
01188 for( long n=0; n<LineSave.nsum; ++n )
01189 {
01190 for( long nEmType=0; nEmType<2; ++nEmType )
01191 {
01192 LineSv[n].SumLine[nEmType+2] += (realnum) LineSv[n].SumLine[nEmType]*cumulative_factor;
01193 }
01194 }
01195
01196 for( long n=0; n<rfield.nflux; ++n)
01197 {
01198
01199 rfield.flux[1][n] += (realnum) rfield.flux[0][n]*cumulative_factor;
01200 rfield.ConEmitReflec[1][n] += (realnum) rfield.ConEmitReflec[0][n]*cumulative_factor;
01201 rfield.ConEmitOut[1][n] += (realnum) rfield.ConEmitOut[0][n]*cumulative_factor;
01202 rfield.ConRefIncid[1][n] += (realnum) rfield.ConRefIncid[0][n]*cumulative_factor;
01203 rfield.flux_total_incident[1][n] += (realnum) rfield.flux_total_incident[0][n]*cumulative_factor;
01204 rfield.reflin[1][n] += (realnum) rfield.reflin[0][n]*cumulative_factor;
01205 rfield.outlin[1][n] += (realnum) rfield.outlin[0][n]*cumulative_factor;
01206 }
01207 }
01208
01209
01210
01211 struc.nzonePreviousIteration = nzone;
01212 for( long i=0; i<struc.nzonePreviousIteration; ++i )
01213 {
01214 struc.depth_last[i] = struc.depth[i];
01215 struc.drad_last[i] = struc.drad[i];
01216 }
01217
01218
01219
01220
01221 for( long i=0; i < rfield.nflux; i++ )
01222 {
01223 {
01224 enum{DEBUG_LOC=false};
01225 if( DEBUG_LOC)
01226 {
01227 fprintf(ioQQQ,"i=%li opac %.2e \n", i,
01228 (double)opac.opacity_abs[i]*radius.drad_x_fillfac/2.*(double)geometry.DirectionalCosin );
01229 }
01230 }
01231 double tau = opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin;
01232 ASSERT( tau > 0. );
01233 fac = sexp( tau );
01234
01235
01236
01237
01238
01239 if( (realnum)(fac/SDIV(rfield.ConInterOut[i]))>SMALLFLOAT && fac > SMALLFLOAT )
01240 {
01241 rfield.ConInterOut[i] /= (realnum)fac;
01242 rfield.outlin[0][i] /= (realnum)fac;
01243 rfield.outlin_noplot[i] /= (realnum)fac;
01244 }
01245 }
01246
01247
01248 radius.StopThickness[iteration-1] = radius.depth;
01249 return;
01250 }