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