cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*IterStart set and save values of many variables at start of iteration */ 00004 /*IterRestart restart iteration */ 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 /* these are the static saved variables, set in IterStart and reset in IterRestart */ 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 /* arguments are atomic number, ionization stage*/ 00062 static realnum xIonFsave[LIMELM+3][LIMELM+1]; 00063 static double HeatSave[LIMELM][LIMELM]; 00064 static realnum supsav[LIMELM][LIMELM], 00065 dndtSave; 00066 /* save values of dr and drNext */ 00067 static double drSave , drNextSave; 00068 00069 /* also places to save them between iterations */ 00070 static long int IonLowSave[LIMELM], 00071 IonHighSave[LIMELM]; 00072 00073 /* these are used to reset the level populations of the h and he model atoms */ 00074 /*static realnum hnsav[LIMELM][LMHLVL+1], */ 00075 static bool lgHNSAV = false; 00076 00077 static realnum ***HOpacRatSav , 00078 H_sum_in_CO_save; 00079 00080 static realnum hmsav[N_H_MOLEC]; 00081 00082 static double ortho_save , para_save; 00083 static double ***hnsav, 00084 edsav; 00085 00086 static realnum xMolFracsSave[LIMELM], 00087 gas_phase_save[LIMELM]; 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 /* allocate two matrices if first time in this core load 00102 * these variables are malloced here because they are file static - 00103 * used to save values at start of first zone, and reset to this value at 00104 * start of new iteration */ 00105 if( !lgHNSAV ) 00106 { 00107 /* set flag so we never do this again */ 00108 lgHNSAV = true; 00109 00110 HOpacRatSav = (realnum ***)MALLOC(sizeof(realnum **)*NISO ); 00111 00112 hnsav = (double ***)MALLOC(sizeof(double **)*NISO ); 00113 00114 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00115 { 00116 00117 HOpacRatSav[ipISO] = (realnum **)MALLOC(sizeof(realnum *)*LIMELM ); 00118 00119 hnsav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ); 00120 00121 /* now do the second dimension */ 00122 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00123 { 00124 HOpacRatSav[ipISO][nelem] = (realnum *)MALLOC(sizeof(realnum)*(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 /* chemistry source and sink terms for ionization ladders */ 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 /* IterStart is called at the start of EVERY iteration */ 00153 if( trace.lgTrace ) 00154 { 00155 fprintf( ioQQQ, " IterStart called.\n" ); 00156 } 00157 00158 /* check if this is the last iteration */ 00159 if( iteration > iterations.itermx ) 00160 { 00161 iterations.lgLastIt = true; 00162 } 00163 else 00164 { 00165 iterations.lgLastIt = false; 00166 } 00167 00168 /* reset some vars dealing with H2 00169 H2_Reset(); */ 00170 00171 /* flag set true in RT_line_one_tauinc if maser ever capped */ 00172 rt.lgMaserCapHit = false; 00173 00174 /* zero out charge transfer heating and cooling fractions */ 00175 atmdat.HCharHeatMax = 0.; 00176 atmdat.HCharCoolMax = 0.; 00177 00178 /* the reason this calculation stops */ 00179 strncpy( StopCalc.chReasonStop, "reason not specified.", 00180 sizeof(StopCalc.chReasonStop) ); 00181 00182 /* zero fractions of He0 destruction due to 23S */ 00183 he.nzone = 0; 00184 he.frac_he0dest_23S = 0.; 00185 he.frac_he0dest_23S_photo = 0.; 00186 00187 /* save expansion rate here */ 00188 dndtSave = dynamics.dDensityDT; 00189 00190 timesc.time_H2_Dest_longest = 0.; 00191 timesc.time_H2_Form_longest = 0.; 00192 /* remains neg if not evaluated */ 00193 timesc.time_H2_Dest_here = -1.; 00194 timesc.time_H2_Form_here = 0.; 00195 00196 for( i=0; i < mole.num_comole_calc; i++ ) 00197 { 00198 timesc.AgeCOMoleDest[i] = 0.; 00199 } 00200 timesc.BigCOMoleForm = 0.; 00201 timesc.TimeH21cm = 0.; 00202 thermal.CoolHeatMax = 0.; 00203 hydro.HCollIonMax = 0.; 00204 00205 atmdat.HIonFracMax = 0.; 00206 00207 /* will save highest n=2p/1s ratio for hydrogen atom*/ 00208 hydro.pop2mx = 0.; 00209 hydro.lgHiPop2 = false; 00210 00211 hydro.nLyaHot = 0; 00212 hydro.TLyaMax = 0.; 00213 00214 /* evaluate the gas and radiation pressure */ 00215 /* this sets values of pressure.PresTotlCurr */ 00216 pressure.PresInteg = 0.; 00217 pressure.pinzon = 0.; 00218 PresTotCurrent(); 00219 00220 dynamics.HeatMax = 0.; 00221 dynamics.CoolMax = 0.; 00222 00223 /* reset these since we now have a valid solution */ 00224 pressure.pbeta = 0.; 00225 pressure.RadBetaMax = 0.; 00226 pressure.lgPradCap = false; 00227 pressure.lgPradDen = false; 00228 /* this flag will say we hit the sonic point */ 00229 pressure.lgSonicPoint = false; 00230 pressure.lgStrongDLimbo = false; 00231 00232 /* define two timescales for equilibrium, Compton and thermal */ 00233 timesc.tcmptn = 1.e0/(rfield.qtot*6.65e-25*dense.eden); 00234 timesc.time_therm_long = 1.5*dense.pden*BOLTZMANN*phycon.te/thermal.ctot; 00235 00236 /* this will be total mass in computed structure */ 00237 dense.xMassTotal = 0.; 00238 00239 for( nelem=0; nelem < LIMELM; nelem++ ) 00240 { 00241 00242 if( dense.lgElmtOn[nelem] ) 00243 { 00244 00245 /* now zero out the ionic fractions */ 00246 for( ion=0; ion < (nelem + 2); ion++ ) 00247 { 00248 xIonFsave[nelem][ion] = dense.xIonDense[nelem][ion]; 00249 } 00250 00251 for( ion=0; ion < LIMELM; ion++ ) 00252 { 00253 HeatSave[nelem][ion] = thermal.heating[ion][nelem]; 00254 } 00255 } 00256 } 00257 00258 /* >>chng 02 jan 28, add ipISO loop */ 00259 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00260 { 00261 /* >>chng 03 apr 11, from 0 to ipISO */ 00262 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00263 { 00264 if( dense.lgElmtOn[nelem] ) 00265 { 00266 for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ ) 00267 { 00268 HOpacRatSav[ipISO][nelem][n] = iso.ConOpacRatio[ipISO][nelem][n]; 00269 hnsav[ipISO][nelem][n] = StatesElem[ipISO][nelem][n].Pop; 00270 } 00271 } 00272 iso.TwoNu_induc_dn_max[ipISO][nelem] = 0.; 00273 } 00274 } 00275 00276 rfield.ipEnergyBremsThin = 0; 00277 rfield.EnergyBremsThin = 0.; 00278 00279 atoms.nNegOI = 0; 00280 for( i=0; i< N_OI_LEVELS; ++i ) 00281 { 00282 atoms.popoi[i] = 0.; 00283 } 00284 00285 /* save molecular fractions and ionization range */ 00286 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00287 { 00288 /* save molecular densities */ 00289 xMolFracsSave[nelem] = dense.xMolecules[nelem]; 00290 gas_phase_save[nelem] = dense.gas_phase[nelem]; 00291 IonLowSave[nelem] = dense.IonLow[nelem]; 00292 IonHighSave[nelem] = dense.IonHigh[nelem]; 00293 } 00294 /*fprintf(ioQQQ,"DEBUG IterStart save set to gas_phase hden %.4e \n", 00295 dense.gas_phase[0]);*/ 00296 00297 edsav = dense.eden; 00298 H_sum_in_CO_save = dense.H_sum_in_CO; 00299 00300 /* save molecular densities */ 00301 for( i=0; i<N_H_MOLEC; i++) 00302 { 00303 hmsav[i] = hmi.Hmolec[i]; 00304 } 00305 ortho_save = h2.ortho_density; 00306 para_save = h2.para_density; 00307 00308 for( i=0; i<mole.num_comole_calc; ++i ) 00309 { 00310 COmole[i]->hevmol_save = COmole[i]->hevmol; 00311 } 00312 00313 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem ) 00314 { 00315 /* these have one more ion than above */ 00316 for( ion=0; ion<nelem+2; ++ion ) 00317 { 00318 long int ion2; 00319 /* zero out the source and sink arrays */ 00320 saveMoleSource[nelem][ion] = mole.source[nelem][ion]; 00321 saveMoleSink[nelem][ion] = mole.sink[nelem][ion]; 00322 /*>>chng 06 jun 27, move from here, where leak occurs, to 00323 * above where xMoleChTrRate first created */ 00324 /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/ 00325 for( ion2=0; ion2<nelem+2; ++ion2 ) 00326 { 00327 SaveMoleChTrRate[nelem][ion][ion2] = mole.xMoleChTrRate[nelem][ion][ion2]; 00328 } 00329 } 00330 } 00331 00332 hmihet_save = hmi.hmihet; 00333 hmitot_save = hmi.hmitot; 00334 00335 h2plus_heat_save = hmi.h2plus_heat; 00336 00337 HeatH2Dish_used_save = hmi.HeatH2Dish_used; 00338 HeatH2Dexc_used_save = hmi.HeatH2Dexc_used; 00339 00340 deriv_HeatH2Dexc_used_save = hmi.deriv_HeatH2Dexc_used; 00341 H2_Solomon_dissoc_rate_used_H2g_save = hmi.H2_Solomon_dissoc_rate_used_H2g; 00342 H2_Solomon_dissoc_rate_used_H2s_save = hmi.H2_Solomon_dissoc_rate_used_H2s; 00343 00344 H2_H2g_to_H2s_rate_used_save = hmi.H2_H2g_to_H2s_rate_used; 00345 00346 H2_photodissoc_used_H2s_save = hmi.H2_photodissoc_used_H2s; 00347 H2_photodissoc_used_H2g_save = hmi.H2_photodissoc_used_H2g; 00348 00349 UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face; 00350 UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_depth; 00351 UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_face; 00352 00353 /* save zone thickness, and next zone thickness */ 00354 drSave = radius.drad; 00355 drNextSave = radius.drNext; 00356 /* remember largest and smallest dr over iteration */ 00357 radius.dr_min_last_iter = BIGFLOAT; 00358 radius.dr_max_last_iter = 0.; 00359 00360 geometry.nprint = iterations.IterPrnt[iteration-1]; 00361 00362 colden.TotMassColl = 0.; 00363 colden.tmas = 0.; 00364 colden.wmas = 0.; 00365 colden.rjnmin = FLT_MAX; 00366 colden.ajmmin = FLT_MAX; 00367 00368 thermal.nUnstable = 0; 00369 thermal.lgUnstable = false; 00370 00371 rfield.extin_mag_B_point = 0.; 00372 rfield.extin_mag_V_point = 0.; 00373 rfield.extin_mag_B_extended = 0.; 00374 rfield.extin_mag_V_extended = 0.; 00375 00376 /* plus 1 is to zero out point where unit integration occurs */ 00377 for( i=0; i < rfield.nupper+1; i++ ) 00378 { 00379 /* diffuse fields continuum */ 00380 /* >>chng 03 nov 08, recover SummedDif */ 00381 rfield.SummedDifSave[i] = rfield.SummedDif[i]; 00382 rfield.ConEmitReflec[i] = 0.; 00383 rfield.ConEmitOut[i] = 0.; 00384 rfield.ConInterOut[i] = 0.; 00385 /* save OTS continuum for next iteration */ 00386 rfield.otssav[i][0] = rfield.otscon[i]; 00387 rfield.otssav[i][1] = rfield.otslin[i]; 00388 rfield.outlin[i] = 0.; 00389 rfield.outlin_noplot[i] = 0.; 00390 rfield.ConOTS_local_OTS_rate[i] = 0.; 00391 rfield.ConOTS_local_photons[i] = 0.; 00392 00393 /* save initial absorption, scattering opacities for next iteration */ 00394 opac.opacity_abs_savzon1[i] = opac.opacity_abs[i]; 00395 opac.opacity_sct_savzon1[i] = opac.opacity_sct[i]; 00396 00397 /* will accumulate optical depths through cloud */ 00398 opac.TauAbsFace[i] = opac.taumin; 00399 opac.TauScatFace[i] = opac.taumin; 00400 /* >>chng 99 dec 04, having exactly 1 zone thickness for first zone caused discontinuity 00401 * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different, 00402 * since tau in is 1e-20, e2 is 0.9999, and so some H ots 00403 * these were not here at all - changed to dr/2 */ 00404 /* attenuation of flux by optical depths IN THIS ZONE 00405 * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command, 00406 * option for illumination of slab at an angle */ 00407 /* >>chng 04 oct 09, from drad to radius.drad_x_fillfac - include fill fac, PvH */ 00408 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin); 00409 00410 /* e(-tau) in inward direction, up to illuminated face */ 00411 opac.ExpmTau[i] = (realnum)opac.ExpZone[i]; 00412 00413 /* e2(tau) in inward direction, up to illuminated face, this is used to get the 00414 * recombination escape probability */ 00415 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i] ); 00416 } 00417 00418 t_fe2ovr_la::Inst().zero_opacity(); 00419 00420 /* this zeros out arrays to hold mean ionization fractions 00421 * later entered by mean 00422 * read out by setlim */ 00423 MeanZero(); 00424 00425 /* zero out the column densities */ 00426 for( i=0; i < NCOLD; i++ ) 00427 { 00428 colden.colden[i] = 0.; 00429 } 00430 colden.He123S = 0.; 00431 colden.H0_ov_Tspin = 0.; 00432 colden.OH_ov_Tspin = 0.; 00433 colden.coldenH2_ov_vel = 0.; 00434 00435 /* upper and lower levels of H0 1s */ 00436 colden.H0_21cm_upper =0; 00437 colden.H0_21cm_lower =0; 00438 00439 h2.ortho_colden = 0.; 00440 h2.para_colden = 0.; 00441 00442 for( i=0; i < mole.num_comole_calc; i++ ) 00443 { 00444 COmole[i]->hevcol = 0.; 00445 } 00446 for( i=0; i < 5; i++ ) 00447 { 00448 colden.C2Pops[i] = 0.; 00449 colden.C2Colden[i] = 0.; 00450 /* pops and column density for SiII atom */ 00451 colden.Si2Pops[i] = 0.; 00452 colden.Si2Colden[i] = 0.; 00453 } 00454 for( i=0; i < 3; i++ ) 00455 { 00456 /* pops and column density for CI atom */ 00457 colden.C1Pops[i] = 0.; 00458 colden.C1Colden[i] = 0.; 00459 /* pops and column density for OI atom */ 00460 colden.O1Pops[i] = 0.; 00461 colden.O1Colden[i] = 0.; 00462 /* pops and column density for CIII atom */ 00463 colden.C3Pops[i] = 0.; 00464 } 00465 for( i=0; i < 4; i++ ) 00466 { 00467 /* pops and column density for CIII atom */ 00468 colden.C3Colden[i] = 0.; 00469 } 00470 00471 /* these are some line of sight emission measures */ 00472 colden.dlnenp = 0.; 00473 colden.dlnenHep = 0.; 00474 colden.dlnenHepp = 0.; 00475 colden.dlnenCp = 0.; 00476 colden.H0_ov_Tspin = 0.; 00477 colden.OH_ov_Tspin = 0.; 00478 00479 /* now zero heavy element molecules */ 00480 molcol("ZERO",ioQQQ); 00481 00482 /* this will be sum of all free free heating over model */ 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.HeatEfficPrimary; 00506 secondaries.savefi = secondaries.SecIon2PrimaryErg; 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 /* these will keep track of the number of convergence failures that occur */ 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 /* abort flag */ 00526 lgAbort = false; 00527 00528 /* zero out averaging routine */ 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 * allocate space for lines arrays 00563 * 00564 ************************************************************************/ 00565 00566 /* there are three types of calls to lines() 00567 * ipass = -1, first call, only count number off lines 00568 * ipass = 0, second pass, save labels and wavelengths 00569 * ipass = 1, integrate intensity*/ 00570 LineSave.ipass = -1; 00571 lines(); 00572 ASSERT( LineSave.nsum > 0 ); 00573 00574 /* in a grid this may not be the first time here, return old memory 00575 * and grab new appropriate for this size, since number of lines to 00576 * be stored can change 00577 * as now coded this memory is freed and reallocated on every iteration. 00578 * this allows some details of line emission to change for last iteration, 00579 * but may not really be necessary */ 00580 if( LineSv!=NULL ) 00581 { 00582 if( trace.lgTrace ) 00583 { 00584 fprintf( ioQQQ, "IterStart has freed LindSv \n" ); 00585 } 00586 free( LineSv ); 00587 LineSv = NULL; 00588 } 00589 00590 /* this is the large main line array */ 00591 LineSv = (LinSv*)MALLOC((unsigned)LineSave.nsum*sizeof( LinSv ) ); 00592 00593 for( i=0; i < LineSave.nsum; i++ ) 00594 { 00595 /* this is array of wavelengths, only defined for some lines, 00596 * negative is sentinel saying whether line should be added into total continuum */ 00597 LineSv[i].wavelength = -1; 00598 } 00599 00600 /* there are three calls to lines() 00601 * first call with ipass = -1 was above, only counted number 00602 * of lines and malloced space. This is second call and will do 00603 * one-time creation of line labels */ 00604 LineSave.ipass = 0; 00605 /* this will count number of lines */ 00606 LineSave.nsum = 0; 00607 lines(); 00608 /* has to be positive */ 00609 ASSERT( LineSave.nsum > 0); 00610 /* in the future calls to lines will result in integrations */ 00611 LineSave.ipass = 1; 00612 00613 if( trace.lgTrace ) 00614 fprintf( ioQQQ, "%7ld lines printed in main line array\n", 00615 LineSave.nsum ); 00616 00617 /* now zero out some variables set by last call to LINES */ 00618 hydro.cintot = 0.; 00619 thermal.totcol = 0.; 00620 rfield.comtot = 0.; 00621 thermal.FreeFreeTotHeat = 0.; 00622 thermal.power = 0.; 00623 thermal.HeatLineMax = 0.; 00624 thermal.GBarMax = 0.; 00625 hyperfine.cooling_max = 0.; 00626 00627 gv.TotalDustHeat = 0.; 00628 gv.dphmax = 0.; 00629 hmi.h2pmax = 0.; 00630 gv.dclmax = 0.; 00631 gv.GrnElecDonateMax = 0.; 00632 gv.GrnElecHoldMax = 0.; 00633 00634 co.codfrc = 0.; 00635 co.codtot = 0.; 00636 00637 hmi.h2dfrc = 0.; 00638 hmi.h2line_cool_frac = 0.; 00639 hmi.h2dtot = 0.; 00640 timesc.sound = 0.; 00641 00642 if( LineSave.lgNormSet ) 00643 { 00644 /* set normalization line index to zero from negative initial value so that 00645 * we pass the assert in cdLine, in order to get the norm line index */ 00646 LineSave.ipNormWavL = 0; 00647 LineSave.ipNormWavL = cdLine( LineSave.chNormLab , LineSave.WavLNorm , &fhe1nx , &ratio ); 00648 if( LineSave.ipNormWavL < 0 ) 00649 { 00650 /* did not find the line if return is negative */ 00651 fprintf( ioQQQ, "PROBLEM could not find the normalisation line.\n"); 00652 fprintf( ioQQQ, 00653 "IterStart could not find a line with a wavelength of %f and a label of %s\n", 00654 LineSave.WavLNorm, LineSave.chNormLab ); 00655 fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n"); 00656 fprintf( ioQQQ, "Sorry.\n"); 00657 cdEXIT(EXIT_FAILURE); 00658 } 00659 } 00660 else 00661 { 00662 /* set normalization line index to zero from negative initial value so that 00663 * we pass the assert in cdLine, in order to get the norm line index */ 00664 LineSave.ipNormWavL = 0; 00665 LineSave.ipNormWavL = cdLine( "TOTL" , 4861. , &fhe1nx , &ratio ); 00666 if( LineSave.ipNormWavL < 0 ) 00667 { 00668 /* did not find the line if return is negative */ 00669 fprintf( ioQQQ, "PROBLEM could not find the default normalisation line.\n"); 00670 fprintf( ioQQQ, 00671 "IterStart could not find a line with a wavelength of 4861 and a label of TOTL\n" ); 00672 fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n"); 00673 fprintf( ioQQQ, "Sorry.\n"); 00674 cdEXIT(EXIT_FAILURE); 00675 } 00676 } 00677 00678 /* set up stop line command on first iteration 00679 * find index for lines and save for future iterations 00680 * StopCalc.nstpl is zero (false) if no stop line commands entered */ 00681 if( iteration == 1 && StopCalc.nstpl ) 00682 { 00683 /* nstpl is number of stop line commands, 0 if none entered */ 00684 for( long int nStopLine=0; nStopLine < StopCalc.nstpl; nStopLine++ ) 00685 { 00686 double relint, absint ; 00687 /* returns array index for line in array stack if we found the line, 00688 * return negative of total number of lines as debugging aid if line not found */ 00689 StopCalc.ipStopLin1[nStopLine] = cdLine( StopCalc.chStopLabel1[nStopLine], 00690 /* wavelength of line in angstroms, not format printed by code */ 00691 StopCalc.StopLineWl1[nStopLine], &relint, &absint ); 00692 00693 if( StopCalc.ipStopLin1[nStopLine]<0 ) 00694 { 00695 fprintf( ioQQQ, 00696 " IterStart could not find first line in STOP LINE command, number %ld with label *%s* and wl %.1f\n", 00697 StopCalc.ipStopLin1[nStopLine] , 00698 StopCalc.chStopLabel1[nStopLine], 00699 StopCalc.StopLineWl1[nStopLine]); 00700 cdEXIT(EXIT_FAILURE); 00701 } 00702 00703 StopCalc.ipStopLin2[nStopLine] = cdLine( StopCalc.chStopLabel2[nStopLine], 00704 /* wavelength of line in angstroms, not format printed by code */ 00705 StopCalc.StopLineWl2[nStopLine], &relint, &absint ); 00706 00707 if( StopCalc.ipStopLin2[nStopLine] < 0 ) 00708 { 00709 fprintf( ioQQQ, 00710 " IterStart could not find second line in STOP LINE command, line number %ld with label *%s* and wl %.1f\n", 00711 StopCalc.ipStopLin1[nStopLine] , 00712 StopCalc.chStopLabel1[nStopLine], 00713 StopCalc.StopLineWl1[nStopLine]); 00714 cdEXIT(EXIT_FAILURE); 00715 } 00716 00717 if( trace.lgTrace ) 00718 { 00719 fprintf( ioQQQ, 00720 " stop line 1 is number %5ld wavelength is %f label is %4.4s\n", 00721 StopCalc.ipStopLin1[nStopLine], 00722 LineSv[StopCalc.ipStopLin1[nStopLine]].wavelength, 00723 LineSv[StopCalc.ipStopLin1[nStopLine]].chALab ); 00724 fprintf( ioQQQ, 00725 " stop line 2 is number %5ld wavelength is %f label is %4.4s\n", 00726 StopCalc.ipStopLin2[nStopLine], 00727 LineSv[StopCalc.ipStopLin2[nStopLine]].wavelength, 00728 LineSv[StopCalc.ipStopLin2[nStopLine]].chALab ); 00729 } 00730 } 00731 } 00732 00733 /* option to only print last iteration */ 00734 if( prt.lgPrtLastIt ) 00735 { 00736 if( iteration == 1 ) 00737 { 00738 called.lgTalk = false; 00739 } 00740 00741 /* initial condition of TALK may be off if AMOEBA used 00742 * sec part is for print last command 00743 * lgTalkForcedOff is set true when cdTalk is called with false 00744 * to turn off printing */ 00745 if( iterations.lgLastIt && !prt.lgPrtStart && !called.lgTalkForcedOff ) 00746 { 00747 called.lgTalk = called.lgTalkSave; 00748 } 00749 } 00750 00751 if( trace.lgTrace && trace.lgOptcBug ) 00752 { 00753 if( iteration > 1 ) 00754 { 00755 fprintf( ioQQQ, " Outer optical depths defined.\n" ); 00756 } 00757 else 00758 { 00759 fprintf( ioQQQ, " Outer optical depths NOT defined.\n" ); 00760 } 00761 00762 /* heavy element lines */ 00763 for( i=0; i < nLevel1; i++ ) 00764 { 00765 fprintf( ioQQQ, "%6f:%8.1e%8.1e%8.1e", 00766 TauLines[i].WLAng, 00767 TauLines[i].Emis->TauIn, 00768 TauLines[i].Emis->TauTot, 00769 TauLines[i].Emis->Pesc ); 00770 } 00771 00772 /*Atomic and molecular lines-Humeshkar Nemala*/ 00773 for( i=0; i <linesAdded2; i++) 00774 { 00775 fprintf( ioQQQ, "%6f:%8.1e%8.1e%8.1e", 00776 atmolEmis[i].tran->WLAng, 00777 atmolEmis[i].TauIn, 00778 atmolEmis[i].TauTot, 00779 atmolEmis[i].Pesc); 00780 } 00781 00782 fprintf( ioQQQ, "\n" ); 00783 } 00784 00785 if( opac.lgCaseB ) 00786 { 00787 if( trace.lgTrace ) 00788 { 00789 fprintf( ioQQQ, " IterStart does not change mid-zone optical depth since CASE B\n" ); 00790 } 00791 } 00792 00793 /* check if induced recombination can be ignored */ 00794 hydro.FracInd = 0.; 00795 hydro.fbul = 0.; 00796 00797 /* remember some things about effects of induced rec on H only 00798 * don't do ground state since SPHERE turns it off */ 00799 for( i=ipH2p; i < (iso.numLevels_max[ipH_LIKE][ipHYDROGEN] - 1); i++ ) 00800 { 00801 if( Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->Aul <= iso.SmallA ) 00802 continue; 00803 00804 ratio = iso.RecomInducRate[ipH_LIKE][ipHYDROGEN][i]*iso.PopLTE[ipH_LIKE][ipHYDROGEN][i]/ 00805 SDIV(iso.RecomInducRate[ipH_LIKE][ipHYDROGEN][i]*iso.PopLTE[ipH_LIKE][ipHYDROGEN][i] + 00806 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][i][ipRecRad]* 00807 iso.RadRecomb[ipH_LIKE][ipHYDROGEN][i][ipRecNetEsc]); 00808 if( ratio > hydro.FracInd ) 00809 { 00810 hydro.FracInd = (realnum)ratio; 00811 hydro.ndclev = i; 00812 } 00813 00814 ratio = Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->pump/ 00815 (Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->pump + 00816 Transitions[ipH_LIKE][ipHYDROGEN][i+1][i].Emis->Aul); 00817 00818 if( ratio > hydro.fbul ) 00819 { 00820 hydro.fbul = (realnum)ratio; 00821 hydro.nbul = i; 00822 } 00823 } 00824 00825 /* this was set in call to lines above */ 00826 ASSERT( LineSave.nsum > 0); 00827 00828 if( trace.lgTrace ) 00829 fprintf( ioQQQ, " IterStart returns.\n" ); 00830 return; 00831 } 00832 00833 /*IterRestart reset many variables at the start of a new iteration 00834 * called by cloudy after calculation is completed, when more iterations 00835 * are needed - the iteration has been incremented when this routine is 00836 * called so iteration == 2 after first iteration, we are starting 00837 * the second iteration */ 00838 void IterRestart(void) 00839 { 00840 long int i, 00841 ion, 00842 ipISO , 00843 n, 00844 nelem; 00845 double SumOTS; 00846 00847 DEBUG_ENTRY( "IterRestart()" ); 00848 00849 dynamics.dDensityDT = dndtSave; 00850 00851 /* this is case where temperature floor has been set, if it was hit 00852 * then we did a constant temperature calculation, and must go back to 00853 * a thermal solution 00854 * test on thermal.lgTemperatureConstantCommandParsed distinguishes 00855 * from temperature floor option, so not reset if constant temperature 00856 * was actually set */ 00857 if( StopCalc.TeFloor > 0. && !thermal.lgTemperatureConstantCommandParsed ) 00858 { 00859 thermal.lgTemperatureConstant = false; 00860 thermal.ConstTemp = 0.; 00861 } 00862 00863 lgAbort = false; 00864 00865 /* reset some parameters needed for magnetic field */ 00866 Magnetic_reinit(); 00867 00868 opac.stimax[0] = 0.; 00869 opac.stimax[1] = 0.; 00870 00871 H2_Reset(); 00872 00873 ca.oldcla = 0.; 00874 00875 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem ) 00876 { 00877 /* these have one more ion than above */ 00878 for( ion=0; ion<nelem+2; ++ion ) 00879 { 00880 long int ion2; 00881 /* zero out the source and sink arrays */ 00882 mole.source[nelem][ion] = saveMoleSource[nelem][ion]; 00883 mole.sink[nelem][ion] = saveMoleSink[nelem][ion]; 00884 /*>>chng 06 jun 27, move from here, where leak occurs, to 00885 * above where xMoleChTrRate first created */ 00886 /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/ 00887 for( ion2=0; ion2<nelem+2; ++ion2 ) 00888 { 00889 mole.xMoleChTrRate[nelem][ion][ion2] = SaveMoleChTrRate[nelem][ion][ion2]; 00890 } 00891 } 00892 } 00893 00894 /* reset molecular abundances */ 00895 for( i=0; i<N_H_MOLEC; i++) 00896 { 00897 hmi.Hmolec[i] = hmsav[i]; 00898 } 00899 hmi.H2_total = hmi.Hmolec[ipMH2g] + hmi.Hmolec[ipMH2s]; 00900 /*fprintf(ioQQQ," IterRestar sets H2 total to %.2e\n",hmi.H2_total );*/ 00901 h2.ortho_density = ortho_save; 00902 h2.para_density = para_save; 00903 00904 /* zero out the column densities */ 00905 for( i=0; i < NCOLD; i++ ) 00906 { 00907 colden.colden[i] = 0.; 00908 } 00909 colden.He123S = 0.; 00910 colden.coldenH2_ov_vel = 0.; 00911 00912 for( i=0; i < mole.num_comole_calc; i++ ) 00913 { 00914 /* column densities */ 00915 COmole[i]->hevcol = 0.; 00916 /* largest fraction of atoms in molecules */ 00917 COmole[i]->xMoleFracMax = 0.; 00918 } 00919 00920 for( i=0; i< mole.num_comole_calc; ++i ) 00921 { 00922 COmole[i]->hevmol = COmole[i]->hevmol_save; 00923 00924 /* check for NaN */ 00925 ASSERT( !isnan( COmole[i]->hevmol ) ); 00926 } 00927 00928 /* close out this iteration if dynamics or time dependence is enabled */ 00929 if( dynamics.lgAdvection ) 00930 DynaEndIter(); 00931 00932 rfield.extin_mag_B_point = 0.; 00933 rfield.extin_mag_V_point = 0.; 00934 rfield.extin_mag_B_extended = 0.; 00935 rfield.extin_mag_V_extended = 0.; 00936 00937 hmi.hmihet = hmihet_save; 00938 hmi.hmitot = hmitot_save; 00939 00940 hmi.BiggestH2 = -1.f; 00941 hmi.h2plus_heat = h2plus_heat_save; 00942 hmi.HeatH2Dish_used = HeatH2Dish_used_save; 00943 hmi.HeatH2Dexc_used = HeatH2Dexc_used_save; 00944 00945 hmi.deriv_HeatH2Dexc_used = (realnum)deriv_HeatH2Dexc_used_save; 00946 hmi.H2_Solomon_dissoc_rate_used_H2g = H2_Solomon_dissoc_rate_used_H2g_save; 00947 hmi.H2_Solomon_dissoc_rate_used_H2s = H2_Solomon_dissoc_rate_used_H2s_save; 00948 hmi.H2_H2g_to_H2s_rate_used = H2_H2g_to_H2s_rate_used_save; 00949 hmi.H2_photodissoc_used_H2s = H2_photodissoc_used_H2s_save; 00950 hmi.H2_photodissoc_used_H2g = H2_photodissoc_used_H2g_save; 00951 00952 hmi.UV_Cont_rel2_Draine_DB96_face = (realnum)UV_Cont_rel2_Draine_DB96_face; 00953 hmi.UV_Cont_rel2_Draine_DB96_depth = (realnum)UV_Cont_rel2_Draine_DB96_depth; 00954 hmi.UV_Cont_rel2_Habing_TH85_face = (realnum)UV_Cont_rel2_Habing_TH85_face; 00955 00956 rfield.ipEnergyBremsThin = 0; 00957 rfield.EnergyBremsThin = 0.; 00958 rfield.lgUSphON = false; 00959 00960 radius.glbdst = 0.; 00961 /* zone thickness, and next zone thickness */ 00962 radius.drad = drSave; 00963 radius.drNext = drNextSave; 00964 00965 /* was min dr ever used? */ 00966 radius.lgDrMinUsed = false; 00967 00968 /* simple estimate of H-beta intensity */ 00969 continuum.cn4861 = continuum.sv4861; 00970 continuum.cn1216 = continuum.sv1216; 00971 00972 /* total luminosity */ 00973 continuum.totlsv = continuum.TotalLumin; 00974 00975 /* set debugger on now if NZONE desired is 0 */ 00976 if( (trace.nznbug == 0 && iteration >= trace.npsbug) && trace.lgTrOvrd ) 00977 { 00978 if( trace.nTrConvg==0 ) 00979 { 00980 trace.lgTrace = true; 00981 } 00982 else 00983 /* trace convergence entered = but with negative flag = make positive, 00984 * abs and not mult by -1 since may trigger more than one time */ 00985 trace.nTrConvg = abs( trace.nTrConvg ); 00986 00987 fprintf( ioQQQ, " IterRestart called.\n" ); 00988 } 00989 else 00990 { 00991 trace.lgTrace = false; 00992 } 00993 00994 /* zero secondary suprathermals variables for first ionization attem */ 00995 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00996 { 00997 for( ion=0; ion<nelem+1; ++ion ) 00998 { 00999 secondaries.csupra[nelem][ion] = supsav[nelem][ion]; 01000 } 01001 } 01002 secondaries.x12tot = secondaries.x12sav; 01003 secondaries.HeatEfficPrimary = secondaries.hetsav; 01004 secondaries.SecIon2PrimaryErg = secondaries.savefi; 01005 for( nelem=0; nelem<LIMELM; ++nelem) 01006 { 01007 for( ion=0; ion<nelem+1; ++ion ) 01008 { 01009 ionbal.CompRecoilHeatRate[nelem][ion] = ionbal.CompRecoilHeatRateSave[nelem][ion]; 01010 ionbal.CompRecoilIonRate[nelem][ion] = ionbal.CompRecoilIonRateSave[nelem][ion]; 01011 } 01012 } 01013 01014 wind.lgVelPos = true; 01015 wind.AccelMax = 0.; 01016 wind.AccelAver = 0.; 01017 wind.acldr = 0.; 01018 wind.windv = wind.windv0; 01019 01020 thermal.nUnstable = 0; 01021 thermal.lgUnstable = false; 01022 01023 pressure.pbeta = 0.; 01024 pressure.RadBetaMax = 0.; 01025 pressure.lgPradCap = false; 01026 pressure.lgPradDen = false; 01027 /* this flag will say we hit the sonic point */ 01028 pressure.lgSonicPoint = false; 01029 pressure.lgStrongDLimbo = false; 01030 pressure.PresInteg = 0.; 01031 pressure.pinzon = 0.; 01032 01033 dense.eden = edsav; 01034 dense.EdenHCorr = edsav; 01035 dense.EdenTrue = edsav; 01036 dense.H_sum_in_CO = H_sum_in_CO_save; 01037 01038 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 01039 { 01040 /* reset molecular densities */ 01041 dense.gas_phase[nelem] = gas_phase_save[nelem]; 01042 dense.xMolecules[nelem] = xMolFracsSave[nelem]; 01043 dense.IonLow[nelem] = IonLowSave[nelem]; 01044 dense.IonHigh[nelem] = IonHighSave[nelem]; 01045 } 01046 /*fprintf(ioQQQ,"DEBUG IterRestart gas_phase set to save hden %.4e\n", 01047 dense.gas_phase[0]);*/ 01048 01049 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 01050 { 01051 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 01052 { 01053 if( dense.lgElmtOn[nelem] ) 01054 { 01055 for( n=ipH1s; n < iso.numLevels_max[ipISO][nelem]; n++ ) 01056 { 01057 iso.ConOpacRatio[ipISO][nelem][n] = HOpacRatSav[ipISO][nelem][n]; 01058 StatesElem[ipISO][nelem][n].Pop = hnsav[ipISO][nelem][n]; 01059 } 01060 } 01061 } 01062 } 01063 01064 if( trace.lgTrace && trace.lgNeBug ) 01065 { 01066 fprintf( ioQQQ, " EDEN set to%12.4e by IterRestart.\n", 01067 dense.eden ); 01068 } 01069 01070 # if 0 01071 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 01072 { 01073 if( dense.lgElmtOn[nelem] ) 01074 { 01075 /* reset total gas phase abundance to the original set */ 01076 dense.gas_phase[nelem] = abund.solar[nelem]*dense.gas_phase[ipHYDROGEN]; 01077 } 01078 else 01079 { 01080 /* >>chng 04 apr 20, set to zero if element is off */ 01081 dense.gas_phase[nelem] = 0.; 01082 } 01083 } 01084 # endif 01085 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 01086 { 01087 for( ion=0; ion < (nelem + 2); ion++ ) 01088 { 01089 dense.xIonDense[nelem][ion] = xIonFsave[nelem][ion]; 01090 } 01091 for( i=0; i < LIMELM; i++ ) 01092 { 01093 thermal.heating[nelem][i] = HeatSave[nelem][i]; 01094 } 01095 } 01096 01097 GrainRestartIter(); 01098 01099 /* continuum was saved in flux_total_incident */ 01100 for( i=0; i < rfield.nupper; i++ ) 01101 { 01102 /* time-constant part of beamed continuum */ 01103 rfield.flux_beam_const[i] = rfield.flux_beam_const_save[i]; 01104 /* continuum flux_time_beam_save has the initial value of the 01105 * time-dependent beamed continuum */ 01106 rfield.flux_beam_time[i] = rfield.flux_time_beam_save[i]* 01107 rfield.time_continuum_scale; 01108 /* the isotropic continuum source is time steady */ 01109 rfield.flux_isotropic[i] = rfield.flux_isotropic_save[i]; 01110 rfield.flux[i] = rfield.flux_isotropic[i] + rfield.flux_beam_time[i] + 01111 rfield.flux_beam_const[i]; 01112 01113 rfield.SummedDif[i] = rfield.SummedDifSave[i]; 01114 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i]; 01115 rfield.trans_coef_zone[i] = 1.f; 01116 rfield.trans_coef_total[i] = 1.f; 01117 01118 rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i]; 01119 rfield.otscon[i] = rfield.otssav[i][0]; 01120 rfield.otslin[i] = rfield.otssav[i][1]; 01121 /* >>chng 01 apr 10, add these two resets */ 01122 /* >>chng 03 may 20, had set to otscon and otslin, which is whatever 01123 * was left from previous iteration - change to zero */ 01124 rfield.otsconNew[i] = 0.; 01125 rfield.otslinNew[i] = 0.; 01126 rfield.ConInterOut[i] = 0.; 01127 rfield.OccNumbBremsCont[i] = 0.; 01128 rfield.OccNumbDiffCont[i] = 0.; 01129 rfield.OccNumbContEmitOut[i] = 0.; 01130 rfield.outlin[i] = 0.; 01131 rfield.outlin_noplot[i] = 0.; 01132 rfield.ConOTS_local_OTS_rate[i] = 0.; 01133 rfield.ConOTS_local_photons[i] = 0.; 01134 01135 opac.opacity_abs[i] = opac.opacity_abs_savzon1[i]; 01136 opac.OldOpacSave[i] = opac.opacity_abs_savzon1[i]; 01137 opac.opacity_sct[i] = opac.opacity_sct_savzon1[i]; 01138 opac.albedo[i] = 01139 opac.opacity_sct[i]/MAX2(1e-30,opac.opacity_sct[i] + opac.opacity_abs[i]); 01140 opac.tmn[i] = 1.; 01141 /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity 01142 * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different, 01143 * since tau in is 1e-20, e2 is 0.9999, and so some H ots 01144 opac.ExpmTau[i] = 1.; 01145 opac.E2TauAbsFace[i] = 1.;*/ 01146 /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity 01147 * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different, 01148 * since tau in is 1e-20, e2 is 0.9999, and so some H ots 01149 * these were not here at all*/ 01150 /* attenuation of flux by optical depths IN THIS ZONE 01151 * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command, 01152 * option for illumination of slab at an angle */ 01153 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad/2.*geometry.DirectionalCosin); 01154 01155 /* e(-tau) in inward direction, up to illuminated face */ 01156 opac.ExpmTau[i] = (realnum)opac.ExpZone[i]; 01157 01158 /* e2(tau) in inward direction, up to illuminated face */ 01159 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]); 01160 rfield.reflin[i] = 0.; 01161 rfield.ConEmitReflec[i] = 0.; 01162 rfield.ConEmitOut[i] = 0.; 01163 rfield.ConRefIncid[i] = 0.; 01164 01165 /* escape in the outward direction 01166 * on second and later iterations define outward E2 */ 01167 if( iteration > 1 ) 01168 { 01169 /* e2 from current position to outer edge of shell */ 01170 realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] ); 01171 opac.E2TauAbsOut[i] = (realnum)e2( tau ); 01172 ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. ); 01173 } 01174 else 01175 opac.E2TauAbsOut[i] = 1.; 01176 01177 } 01178 01179 /* update continuum */ 01180 RT_OTS_Update(&SumOTS , 0.); 01181 01182 thermal.FreeFreeTotHeat = 0.; 01183 01184 if( called.lgTalk ) 01185 { 01186 fprintf( ioQQQ, "\f\n Start Iteration Number%2ld %75.75s\n\n\n", 01187 iteration, input.chTitle ); 01188 } 01189 01190 /* reset variable to do with FeII atom, in FeIILevelPops */ 01191 FeIIReset(); 01192 return; 01193 } 01194 01195 /* do some work with ending the iteration */ 01196 void IterEnd(void) 01197 { 01198 long i; 01199 double fac, 01200 tau; 01201 01202 DEBUG_ENTRY( "IterEnd()" ); 01203 01204 /* give indication of geometry */ 01205 fac = radius.depth/radius.rinner; 01206 if( fac < 0.1 ) 01207 { 01208 geometry.lgGeoPP = true; 01209 } 01210 else 01211 { 01212 geometry.lgGeoPP = false; 01213 } 01214 01215 /* >>chng 06 apr 03, add last radius and drad */ 01216 for( i=0; i<struc.nzone; ++i ) 01217 { 01218 struc.depth_last[i] = struc.depth[i]; 01219 struc.drad_last[i] = struc.drad[i]; 01220 } 01221 struc.nzone_last = struc.nzone; 01222 01223 /* all continue were attenuated in last zone in radius_increment to represent the intensity 01224 * in the middle of the next zone - this was too much since we now want 01225 * intensity emergent from outer edge of last zone */ 01226 for( i=0; i < rfield.nflux; i++ ) 01227 { 01228 { 01229 enum{DEBUG_LOC=false}; 01230 if( DEBUG_LOC) 01231 { 01232 fprintf(ioQQQ,"i=%li opac %.2e \n", i, 01233 (double)opac.opacity_abs[i]*radius.drad_x_fillfac/2.*(double)geometry.DirectionalCosin ); 01234 } 01235 } 01236 tau = opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin; 01237 ASSERT( tau > 0. ); 01238 fac = sexp( tau ); 01239 01240 /* >>chng 02 dec 14, add first test to see whether product of ratio is within 01241 * range of a float - ConInterOut can be finite and fac just above a small float, 01242 * so ratio exceeds largest size of a float */ 01243 /*if( fac > SMALLFLOAT )*/ 01244 if( (realnum)(fac/SDIV(rfield.ConInterOut[i]))>SMALLFLOAT && fac > SMALLFLOAT ) 01245 { 01246 rfield.ConInterOut[i] /= (realnum)fac; 01247 rfield.outlin[i] /= (realnum)fac; 01248 rfield.outlin_noplot[i] /= (realnum)fac; 01249 } 01250 } 01251 01252 /* remember thickness of previous iteration */ 01253 radius.router[iteration-1] = radius.depth; 01254 return; 01255 }