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