/*StartIter set and save values of many variables at start of iteration */ /*RestartIter restart iteration */ #include #include #include #include "cddefines.h" #include "physconst.h" #include "nhe1lvl.h" #include "nhe3lvl.h" #include "taulines.h" #include "hydrogenic.h" #include "facexp.h" #include "sexp.h" #include "angllum.h" #include "e2.h" #include "elmton.h" #include "grainvar.h" #include "pop371.h" #include "pressure.h" #include "linesave.h" #include "stopcalc.h" #include "fe2neg.h" #include "converge.h" #include "lyaext.h" #include "lamase.h" #include "mean.h" #include "densty.h" #include "chmax.h" #include "hemase.h" #include "hhe2phtots.h" #include "comneg.h" #include "hevmolec.h" #include "lotsco.h" #include "heots.h" #include "negoi.h" #include "tff.h" #include "wind.h" #include "colmas.h" #include "he3pht.h" #include "nustbl.h" #include "opac.h" #include "qs.h" #include "timesc.h" #include "frcind.h" #include "max2pht.h" #include "phe1lv.h" #include "he1nxt.h" #include "tehigh.h" #include "he3tau.h" #include "heopfr.h" #include "ophf.h" #include "hphoto.h" #include "dndt.h" #include "he1tau.h" #include "dalerr.h" #include "trace.h" #include "zonecnt.h" #include "itercnt.h" #include "taumin.h" #include "tcool.h" #include "heat.h" #include "dhla.h" #include "colden.h" #include "norm.h" #include "abrems.h" #include "deptau.h" #include "secondaries.h" #include "hmi.h" #include "compton.h" #include "radius.h" #include "h.h" #include "he1bn.h" #include "he3bn.h" #include "he.h" #include "phycon.h" #include "tlast.h" #include "plast.h" #include "touton.h" #include "caseb.h" #include "doppvel.h" #include "called.h" #include "the1lo.h" #include "dphoht.h" #include "codish.h" #include "h2dish.h" #include "intalk.h" #include "pstart.h" #include "he1stat.h" #include "ionfracs.h" #include "heating.h" #include "heatlmax.h" #include "dcala.h" #include "extramax.h" #include "elecden.h" #include "sumoutinc.h" #include "rfield.h" #include "ionrange.h" #include "charexc.h" #include "cap4.h" #include "lines.h" #include "aver.h" #include "molcol.h" #include "globul.h" #include "h2max.h" #include "cmfrac.h" #include "jmin.h" #include "smbeta.h" #include "totlum.h" #include "hmihet.h" #include "hlmax.h" #include "input.h" #include "drminu.h" #include "negdrg.h" #include "uspher.h" #include "dtmase.h" #include "exptau.h" #include "e2tau.h" #include "solar.h" #include "flxsav.h" #include "trovrd.h" #include "pca2ex.h" #include "stimax.h" #include "sumcontinuum.h" #include "startenditer.h" /* these are the static saved variables, set in StartIter and reset in RestartIter */ /* args are atomic number, ionization stage*/ static float xIonFsave[LIMELM + 2][LIMELM]; static double HeatSave[LIMELM][LIMELM]; static float he1sav, he2sav, he3sav, heisav; /* these are used to reset the level popls of the h and he model atoms */ /*static float hnsav[LIMELM][LMHLVL+1], */ static int lgHNSAV = FALSE; static float **hnsav, he1nsav[NHE1LVL + 1]; /*save opacity ratios for hydrogen and heliumm*/ /*float HOpacRatSav[LIMELM][LMHLVL+1], */ float **HOpacRatSav, ophe1s[NHE1LVL]; static double heb1sv[10], heb3sv[NHE3LVL]; static float h2sav, h2psav, hmisav, h3psav; static float hisv, hiisv, hsv, edsav, EDOthSav, EDMSav; void StartIter(void) { char chCAPS[5]; int lgHit; long int i, ii, ion, ipZ, j, n, nd; double dtau, fhe1nx, ratio; # ifdef DEBUG_FUN fputs( "<+>StartIter()\n", debug_fp ); # endif /* allocate two matrices if first time in this coreload */ if( !lgHNSAV ) { lgHNSAV = TRUE; if( (HOpacRatSav = (float **)malloc(sizeof(float *)*LIMELM ))==NULL ) { fprintf( ioQQQ, " could not malloc HOpacRatSav arrays in 1D\n" ); puts( "[Stop in StartIter]" ); exit(1); } if( (hnsav = (float **)malloc(sizeof(float *)*LIMELM ))==NULL ) { fprintf( ioQQQ, " could not malloc hnsav arrays in 1D\n" ); puts( "[Stop in StartIter]" ); exit(1); } /* now do the second dimension */ for( i=0; i IterCnt.itermx ) { IterCnt.lgLastIt = TRUE; } else { IterCnt.lgLastIt = FALSE; } hemase.lgHeMase = FALSE; SumOutInc.SumOutCon = 0.; SumOutInc.SumIncCon = 0.; /* two photon absorption deep in cloud */ hhe2phtots.H12phtOTS = 0.; hhe2phtots.He12phtOTS = 0.; hhe2phtots.He22phtOTS = 0.; /* flag set true in LineTauInc if maser ever capped */ dtMase.lgMaserCapHit = FALSE; /* zero out charge transfer heating and cooling fractions */ CharExc.HCharHeatMax = 0.; CharExc.HCharCoolMax = 0.; /* save exansion rate here */ dndt.dndtSave = dndt.dDensityDT; timesc.AgeH2MoleDest = 0.; timesc.AgeCOMoleDest = 0.; timesc.BigH2MoleForm = 0.; timesc.BigCOMoleForm = 0.; CHMax.CoolHeatMax = 0.; /* will save highest n=2p/1s ratio for hydrogen atom*/ hydro.pop2mx = 0.; hydro.lgHiPop2 = FALSE; LyaExT.nLyaHot = 0; LyaExT.TLyaMax = 0.; Max2Pht.xMax2Pht = 0.; LaMase.xLaMase = 0.; Fe2Neg.nFe2Neg = 0; /* TLAST is the temperature of the last zone, and is normally * evaluated in ZoneDone, at the end of a zone calculation * it has been undefined up til now, and for the first zone is nonsense */ tlastCom.tlast = phycon.te; /* evaluate the gas and radiation pressure */ pressure.PressureInit = (float)TotalPressure(); pressure.pinit = pressure.PressureInit; pressure.PresInteg = 0.; pressure.pinzon = 0.; /* reset these since we now have a valid solution */ pressure.pbeta = 0.; pressure.RadBetaMax = 0.; pressure.lgPradCap = FALSE; pressure.lgPradDen = FALSE; lotsco.lgLotsCO = FALSE; comneg.lgComNeg = FALSE; /* define two timescales for equilibrium, Compton and thermal */ timesc.tcmptn = 1.e0/(qs.qtot*6.65e-25*phycon.eden); timesc.ttherm = 1.5*densty.pden*BOLTZMANN*phycon.te/tcool.ctot; for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { for( n=IP1S; n <= nhlevel; n++ ) { HOpacRatSav[ipZ][n] = ophf.HOpacRatio[ipZ][n]; hnsav[ipZ][n] = hydro.hn[ipZ][n]; HPhoto.cinduc[ipZ][n] = 0.; } /* now zero out the ionic fractions */ for( ion=1; ion <= (ipZ + 2); ion++ ) { xIonFsave[ion][ipZ] = xIonFracs[ion][ipZ]; } for( i=0; i < LIMELM; i++ ) { HeatSave[i][ipZ] = HeatingCom.heating[i][ipZ]; } } } tffCom.ntff = 1; tffCom.tff = 0.; negoi.nNegOI = 0; /* save ionization bal for second row elements */ for( i=0; i < LIMELM; i++ ) { IonRange.IonLowSave[i] = IonRange.IonLow[i]; IonRange.IonHighSave[i] = IonRange.IonHigh[i]; } hsv = (float)phycon.hden; edsav = (float)phycon.eden; EDOthSav = ElecDen.EdenNotH; EDMSav = ElecDen.EdenMet; /* save molecular fractions */ for( i=0; i < NHEVML; i++ ) { hevmolec.HevMolSav[i] = hevmolec.hevmol[i]; } h2sav = hmi.htwo; h2psav = hmi.h2plus; h3psav = hmi.h3plus; hmisav = hmi.hminus; /* save zone thicknes, and next zone thickness */ radius.drSave = (double)radius.drad; radius.drnsv = radius.drNext; /* these are opacity factors needed for recombination coef */ for( n=0; n < (NHE1LVL + 1); n++ ) { he1nsav[n] = phe1lv.he1n[n]; heb1sv[n] = he1bnCOM.he1bn[n]; } for( n=0; n < NHE1LVL; n++ ) { ophe1s[n] = heopfr.ophe1f[n]; } ZoneCnt.nprint = IterCnt.IterPrnt[IterCnt.iter-1]; hisv = h.hi; hiisv = h.hii; colmas.TotMassColl = 0.; colmas.tmas = 0.; colmas.wmas = 0.; the1loCom.nhe1lo = 0; nustbl.nUnstable = 0; nustbl.lgUnstable = FALSE; heisav = he.heiii; he2sav = he.heii; he1sav = he.hei; he.hei = he.hei1; he3sav = he.hei3; for( i=0; i < NHE3LVL; i++ ) { heb3sv[i] = he3bnCom.he3bn[i]; } he3pht.he3photo = 0.; he3pht.he3lya = 0.; for( i=0; i < NCELL; i++ ) { /* diffuse fields continuum */ rfield.ConOutNoInter[i] = 0.; rfield.ConRefNoInter[i] = 0.; /* save OTS continuum for next iteration */ rfield.otssav[0][i] = rfield.otscon[i]; rfield.otssav[1][i] = rfield.otslin[i]; rfield.outcon[i] = 0.; rfield.outlin[i] = 0.; /* save opacities for next iteration */ opac.opacsv[i] = opac.opac[i]; opac.sctsav[i] = opac.scatop[i]; /* will accumulate optical depths through cloud */ deptau.depabs[i] = tauminCom.taumin; deptau.depsct[i] = tauminCom.taumin; /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity * for heating in very high T models in map.in. zone 1 and 2 were 20% different, * since tau in is 1e-20, e2 is 0.9999, and so some H ots * these were not here at all*/ /* attenuation of flux by optical depths IN THIS ZONE * AngleIllum is 1/COS(theta), is usually 1, reset with illuminate command, * option for illumination of slab at an angle */ facexp.ExpZone[i] = sexp(opac.opac[i]*radius.drad/2.*Angllum.AngleIllum); /* e(-tau) in inward direction, up to illuminated face */ exptau.ExpmTau[i] = (float)facexp.ExpZone[i]; /* e2(tau) in inward direction, up to illuminated face, this is used to get the * recombination escape probability */ e2tau.e2TauAbs[i] = (float)e2(deptau.depabs[i],exptau.ExpmTau[i]); } /* this zeroes out arrays to hold mean ionization fractions * later enered by mean * read out by setlim */ MeanZero(); /* zero out the column densities */ for( i=0; i < NCOLD; i++ ) { coldenCom.colden[i] = 0.; } /* now zero heavy element molecules */ molcol("ZERO"); /* this will be sum of all free free heating over model */ hydro.FreeFreeTotHeat = 0.; tehigh.thist = 0.; tehigh.tlowst = 1e20f; wind.AccelAver = 0.; wind.acldr = 0.; dalerr.ilt = 0; dalerr.iltln = 0; dalerr.ilthn = 0; dalerr.ihthn = 0; dalerr.ifail = 0; Secondaries.supsav = Secondaries.csupra; Secondaries.x12sav = Secondaries.x12tot; Secondaries.hetsav = Secondaries.heatef; Secondaries.savefi = Secondaries.efionz; /* these will keep track of the number of convergence failures that occur */ conv.nTeFail = 0; conv.nPreFail = 0; conv.failmx = 0.; conv.nIonFail = 0; conv.nNeFail = 0; conv.nGrainFail = 0; dhla.dstpow = 0.; dhla.dstotl = 0.; dhla.dstsum = 0.; dhla.DustHeatLya = 0.; dhla.dincon = 0.; dhla.HDustDif = 0.; GrainVar.thtot = 0.; abrems.dlnenp = 0.; GrainVar.reftot = 0.; GrainVar.avhden = 0.; /* zero out averaging routine */ aver("zero",0.,0.," "); for( nd=0; nd < NDUST; nd++ ) { /* >>chng 97 july 5, save and reset this * save grain potential */ GrainVar.dstpotsav[nd] = GrainVar.dstpot[nd]; GrainVar.avdust[nd] = 0.; GrainVar.avdpot[nd] = 0.; GrainVar.avdft[nd] = 0.; GrainVar.TeGrainMax[nd] = -1.f; } Compton.comtot = 0.; codish.codfrc = 0.; codish.codtot = 0.; h2dish.h2dfrc = 0.; h2dish.h2dtot = 0.; dphoht.TotalDustHeat = 0.; dphoht.dphmax = 0.; dphoht.h2pmax = 0.; dphoht.dclmax = 0.; heatlmax.HeatLineMax = 0.; ExtraMax.GBarMax = 0.; HPhoto.cintot = 0.; ZoneCnt.lgZoneTrp = FALSE; /************************************************************************ * * allocate space for lines arrays * ************************************************************************/ /* set up line array */ /* first count number of lines */ LineSave.ipass = -1; lines(); /* in a grid this may not be the first time here, return old memory * and grab new appropriate for this size, since number of lines to * be stored can change * as now coded this memory is freed and reallocated on every iteration. * this allows some details of line emission to change for last iteration, * but may not really be necessary */ if( LineDSv!=NULL ) { if( trace.lgTrace ) { fprintf( ioQQQ, "StartIter has freed LindDSv \n" ); } free(LineDSv ); LineDSv = NULL; } if( LineSv!=NULL ) { if( trace.lgTrace ) { fprintf( ioQQQ, "StartIter has freed LindSv \n" ); } free(LineSv ); LineSv = NULL; } /* only allocate space if this sum is positive - will be zero if no grains present */ if( LineSave.ndsum > 0 ) { if( (LineDSv=(LinDstSv *)malloc((unsigned)LineSave.ndsum*sizeof( LinDstSv ) ) ) == NULL ) { fprintf( ioQQQ, "StartIter could not allocate space for LineSave.ndsum array, needed %li\n", LineSave.ndsum ); puts( "[Stop in StartIter]" ); exit(1); } } /* this is the large main line array */ if( (LineSv=(LinSv*)malloc((unsigned)LineSave.nsum*sizeof( LinSv ) ) ) == NULL ) { fprintf( ioQQQ, "StartIter could not allocate space for LineSave.nsum array, needed %li\n", LineSave.nsum ); puts( "[Stop in StartIter]" ); exit(1); } /* this is array of line energies, only defined for some lines */ for( i=0; i < LineSave.nsum; i++ ) { LineSv[i].xLineEnergy = 0.; } LineSave.nsum = 0; /* second call to set up line labels */ LineSave.ipass = 0; lines(); /* in the future calls to lines will result in integrations */ LineSave.ipass = 1; if( trace.lgTrace ) { fprintf( ioQQQ, "%7ld lines printed in main line array\n", LineSave.nsum ); } /* now zero out some variables set by last call to LINES */ HPhoto.cintot = 0.; tcool.totcol = 0.; tcool.GrossHeat = 0.; Compton.comtot = 0.; heat.power = 0.; dphoht.TotalDustHeat = 0.; dphoht.dphmax = 0.; dphoht.h2pmax = 0.; dphoht.dclmax = 0.; heatlmax.HeatLineMax = 0.; ExtraMax.GBarMax = 0.; codish.codfrc = 0.; codish.codtot = 0.; h2dish.h2dfrc = 0.; h2dish.h2dtot = 0.; timesc.sound = 0.; /* find which line to normalise array to (default is H BETA) * NormWL=1 then use first line in list (H-BETA) * NormWL<0 then find line with wavelength=-NormWL */ if( norm.ipNormWavL < 0 ) { lgHit = FALSE; /* repair norm line wavelength */ norm.ipNormWavL = -norm.ipNormWavL; i = 0; /* this must be reset or it will crash */ j = -LONG_MAX; while( !lgHit && i < LineSave.nsum ) { if( LineSv[i].lin == norm.ipNormWavL ) { /* no label was set, so any will do */ if( strcmp( norm.chNormLab," ") == 0 ) { j = i; lgHit = TRUE; } else { /* check whether labels match */ cap4(chCAPS , (char*)LineSv[i].chALab); if( strcmp( norm.chNormLab,chCAPS) == 0 ) { /* when chNormLab set, only count as hit when label matches too */ j = i; lgHit = TRUE; } } } ++i; } if( !lgHit ) { /* did not find the line */ fprintf( ioQQQ, " couldnt find%5ld lab=%4.4s\n", norm.ipNormWavL, norm.chNormLab ); puts( "[Stop in StartIter]" ); exit(1); } /* following confirms that index within bounds, mostly for lint */ assert( j>0 ); assert( j>chng 96 dec 12, got rid of if, change to do while */ while( !lgHit && i < LineSave.nsum ) { if( LineSv[i-1].lin == StopCalc.LineStopWl[ii] ) { StopCalc.ipStopLin1[ii] = i; lgHit = TRUE; } i += 1; } if( !lgHit ) { fprintf( ioQQQ, " STARTER2 couldnt find line%5ld\n", StopCalc.ipStopLin1[ii] ); puts( "[Stop in StartIter]" ); exit(1); } if( trace.lgTrace ) { fprintf( ioQQQ, " stop line 1 is number%5ld wavelength is%9ld label is %4.4s\n", StopCalc.ipStopLin1[ii], LineSv[StopCalc.ipStopLin1[ii]-1].lin, LineSv[StopCalc.ipStopLin1[ii]-1].chALab ); } if( StopCalc.ipStopLin2[ii] != 1 ) { i = 1; lgHit = FALSE; while( !lgHit && i < LineSave.nsum ) { /* >>chng 96 dec 12, more to do while logic * do i=1,nsum */ if( LineSv[i-1].lin == StopCalc.ipStopLin2[ii] ) { StopCalc.ipStopLin2[ii] = i; lgHit = TRUE; } i += 1; } if( !lgHit ) { fprintf( ioQQQ, " STARTER3 couldnt find line%5ld\n", StopCalc.ipStopLin2[ii] ); puts( "[Stop in StartIter]" ); exit(1); } if( trace.lgTrace ) { fprintf( ioQQQ, " stop line 2 is number%5ld wavelength is%9ld label is %4.4s\n", StopCalc.ipStopLin2[ii], LineSv[StopCalc.ipStopLin2[ii]-1].lin, LineSv[StopCalc.ipStopLin2[ii]-1].chALab ); } } } } /* option to only print last iteration */ if( plast.lgPrtLastIt ) { if( IterCnt.iter == 1 ) { called.lgTalk = FALSE; } /* initial condition of TALK may be off if AMOEBA used * sec part is for print last command * lgTalkForce is set true when cdTalk is called with false * to turn off printing */ if( IterCnt.lgLastIt && !pstart.lgPrtStart && !called.lgTalkForce ) { called.lgTalk = intalk.lgInTalk; } } if( trace.lgTrace && trace.lgOptcBug ) { if( touton.lgTauOutOn ) { fprintf( ioQQQ, " Outer optical depths defined.\n" ); } else { fprintf( ioQQQ, " Outer optical depths NOT defined.\n" ); } fprintf( ioQQQ, " Helium \n" ); for( i=1; i <= NHE3TAU; i++ ) { fprintf( ioQQQ, "%6f:%8.1e%8.1e%8.1e", he3tau[i-1].WLAng, he3tau[i-1].TauIn, he3tau[i-1].TauTot, he3tau[i-1].Pesc ); } fprintf( ioQQQ, "\n" ); /* heavy element lines */ for( i=1; i <= nLevel1; i++ ) { fprintf( ioQQQ, "%6f:%8.1e%8.1e%8.1e", TauLines[i-1].WLAng, TauLines[i-1].TauIn, TauLines[i-1].TauTot, TauLines[i-1].Pesc ); } fprintf( ioQQQ, "\n" ); } if( caseb.lgCaseb ) { if( trace.lgTrace ) { fprintf( ioQQQ, " StartIter does not change mid-zone optical depth since CASE B\n" ); } } else { fhe1nx = he.heii*radius.dReff/DoppVel.doppler[1]*0.52; for( j=1; j < NHE1LVL; j++ ) { dtau = (phe1lv.he1n[0] - phe1lv.he1n[j]*he1statCOM.he1stat[0]/ he1statCOM.he1stat[j])*he1tauCOM.he1opc[0][j]; he1nxtCOM.he1nxt[0][j] = (float)MAX2(he1tauCOM.he1tau[0][j], dtau*fhe1nx); } /* shit */ for( i=1; i < (NHE1LVL - 1); i++ ) { for( j=i + 1; j < NHE1LVL; j++ ) { dtau = (phe1lv.he1n[i] - phe1lv.he1n[j]*he1statCOM.he1stat[i]/ he1statCOM.he1stat[j])*he1tauCOM.he1opc[i][j]; he1nxtCOM.he1nxt[i][j] = (float)MAX2(tauminCom.taumin,dtau* fhe1nx); } } if( trace.lgTrace ) { if( trace.lgHe1Bug ) { fprintf( ioQQQ, " StartIter estimates optical depths for center of first zone; HE1NXT array follows:\n" ); for( i=1; i < NHE1LVL; i++ ) { fprintf( ioQQQ, " %1ld", i ); for( j=0; j < 7; j++ ) { fprintf( ioQQQ, "%10.2e", he1nxtCOM.he1nxt[j][i] ); } fprintf( ioQQQ, "\n" ); } } } } /* check if induced recombination can be ignored */ frcind.FracInd = 0.; frcind.fbul = 0.; /* remember some things about effects of induced rec on H onlye * don't do ground state since SPHERE turns it off */ for( i=IP2P; i <= (nhlevel - 1); i++ ) { ratio = HPhoto.rinduc[0][i]/(HPhoto.rinduc[0][i] + hydro.hrec[0][i][ipRecRad]* hydro.hrec[0][i][ipRecNetEsc]); if( ratio > frcind.FracInd ) { frcind.FracInd = (float)ratio; frcind.ndclev = i; } ratio = HydroLines[0][i+1][i].pump/ (HydroLines[0][i+1][i].pump + HydroLines[0][i+1][i].Aul); if( ratio > frcind.fbul ) { frcind.fbul = (float)ratio; frcind.nbul = i; } } if( trace.lgTrace ) { fprintf( ioQQQ, " StartIter returns.\n" ); } # ifdef DEBUG_FUN fputs( " <->StartIter()\n", debug_fp ); # endif return; } /*===============================================================*/ /*reset reset many variables at the start of a new iteration */ void RestartIter(void) { long int i, ion, ipZ, n, ipOTSchange, nd, nelem; double SumOTS , RatioOTSInc; # ifdef DEBUG_FUN fputs( "<+>RestartIter()\n", debug_fp ); # endif SumOutInc.SumOutCon = 0.; SumOutInc.SumIncCon = 0.; dndt.dDensityDT = dndt.dndtSave; h2max.BiggestH2 = -1.f; called.lgBusted = FALSE; jmin.rjnmin = FLT_MAX; jmin.ajmmin = FLT_MAX; LaMase.xLaMase = 0.; stimaxCom.stimax[0] = 0.; stimaxCom.stimax[1] = 0.; dcala.oldcla = 0.; cmfrac.CarMolFrac = -1.f; /* reset molecular fractions */ for( i=0; i < NHEVML; i++ ) { hevmolec.hevmol[i] = hevmolec.HevMolSav[i]; } hmi.htwo = h2sav; hmi.h2plus = h2psav; hmi.h3plus = h3psav; hmi.hminus = hmisav; tffCom.ntff = 1; tffCom.tff = 0.; globul.glbdst = 0.; /* zone thickness, and next zone thickness */ radius.drad = radius.drSave; radius.drNext = radius.drnsv; uspher.lgUSphON = FALSE; drminu.lgDrMinUsed = FALSE; negdrg.lgNegGrnDrg = FALSE; lotsco.lgLotsCO = FALSE; /* simple estimate of H-beta intensity */ smbeta.SimHBeta = (float)(qs.qhtot*4.75e-13); smbeta.cn4861 = smbeta.sv4861; smbeta.cn1216 = smbeta.sv1216; /* total luminosity */ totlum.totlsv = totlum.TotalLumin; /* set debugger on now if NZONE desired is 0 */ if( (trace.nznbug == 0 && IterCnt.iter >= trace.npsbug) && trovrd.lgTrOvrd ) { trace.lgTrace = TRUE; fprintf( ioQQQ, " RestartIter called.\n" ); } else { trace.lgTrace = FALSE; } /* zero secondary suptrathermals variables for first ionization attem */ Secondaries.csupra = Secondaries.supsav; Secondaries.x12tot = Secondaries.x12sav; Secondaries.heatef = Secondaries.hetsav; Secondaries.efionz = Secondaries.savefi; wind.lgVelPos = TRUE; wind.AccelMax = 0.; wind.AccelAver = 0.; wind.acldr = 0.; wind.windv = wind.windv0; nustbl.nUnstable = 0; nustbl.lgUnstable = FALSE; pressure.pbeta = 0.; pressure.RadBetaMax = 0.; pressure.lgPradCap = FALSE; pressure.lgPradDen = FALSE; pca2ex.popca2ex = 0.; /* these are opacity factors needed for recombination coef */ for( n=0; n < NHE1LVL; n++ ) { heopfr.ophe1f[n] = ophe1s[n]; } for( n=0; n < (NHE1LVL + 1); n++ ) { phe1lv.he1n[n] = he1nsav[n]; } heots.esc626 = 0.; heots.esc584 = 0.; heots.esc660 = 0.; heots.esc911 = 0.; heots.esc228 = 0.; for( i=1; i <= nhlevel; i++ ) { /* hlmax[1] is 1s */ hlmaxCom.hlmax[i] = 0.; } phycon.hden = hsv; phycon.eden = edsav; ElecDen.EdenHCorr = (float)phycon.eden; ElecDen.EdenNotH = EDOthSav; ElecDen.EdenMet = EDMSav; ElecDen.EdenTrue = edsav; for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { /* turn H matrix back on */ for( n=IP1S; n <= nhlevel; n++ ) { ophf.HOpacRatio[ipZ][n] = HOpacRatSav[ipZ][n]; HPhoto.cinduc[ipZ][n] = 0.; hydro.hn[ipZ][n] = hnsav[ipZ][n]; } } } /* NHE1LVL is 9, with 2s at 10, this loop is through 10 */ for( n=0; n < NHE1LVL+1; n++ ) { he1bnCOM.he1bn[n] = heb1sv[n]; } for( n=0; n < NHE3LVL; n++ ) { he3bnCom.he3bn[n] = heb3sv[n]; } if( trace.lgTrace && trace.lgNeBug ) { fprintf( ioQQQ, " EDEN set to%12.4e by RestartIter.\n", phycon.eden ); } for( nelem=0; nelem < LIMELM; nelem++ ) { xIonFracs[0][nelem] = (float)(solarCom.solar[nelem]*phycon.hden); for( ion=1; ion <= (nelem + 2); ion++ ) { xIonFracs[ion][nelem] = xIonFsave[ion][nelem]; } for( i=0; i < LIMELM; i++ ) { HeatingCom.heating[i][nelem] = HeatSave[i][nelem]; } } h.hi = hisv; h.hii = hiisv; he.ahe = (float)(solarCom.solar[1]*phycon.hden); he.heiii = heisav; he.heii = he2sav; he.hei = he1sav; he.hei1 = he.hei; he.hei3 = he3sav; hmihetCom.hmihet = 0.; hmihetCom.hmitot = 0.; dhla.DustHeatLya = 0.; dhla.dhcol = 0.; dhla.dincon = 0.; dhla.HDustDif = 0.; dhla.dstpow = 0.; dhla.dstsum = 0.; dhla.dstotl = 0.; for( nd=0; nd < NDUST; nd++ ) { /* >>chng 97 july 5, reset grain potential * reset grain to pential to initial value from previous iteration */ GrainVar.dstpot[nd] = GrainVar.dstpotsav[nd]; GrainVar.avdust[nd] = 0.; GrainVar.avdpot[nd] = 0.; GrainVar.avdft[nd] = 0.; GrainVar.TeGrainMax[nd] = -1.f; } /* fix number of stages of ionization */ for( i=0; i < LIMELM; i++ ) { IonRange.IonLow[i] = IonRange.IonLowSave[i]; IonRange.IonHigh[i] = IonRange.IonHighSave[i]; } /* continuum was saved in FluxSave */ for( i=0; i < NCELL; i++ ) { rfield.flux[i] = flxsav.FluxSave[i]; rfield.occnum[i] = rfield.flux[i]*rfield.convoc[i]; rfield.otscon[i] = rfield.otssav[0][i]; rfield.otslin[i] = rfield.otssav[1][i]; rfield.outcon[i] = 0.; rfield.outlin[i] = 0.; opac.opac[i] = opac.opacsv[i]; opac.OldOpacSave[i] = opac.opacsv[i]; opac.scatop[i] = opac.sctsav[i]; opac.albedo[i] = opac.scatop[i]/MAX2(1e-30,opac.scatop[i] + opac.opac[i]); opac.tmn[i] = 1.; /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity * for heating in very high T models in map.in. zone 1 and 2 were 20% different, * since tau in is 1e-20, e2 is 0.9999, and so some H ots exptau.ExpmTau[i] = 1.; e2tau.e2TauAbs[i] = 1.;*/ /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity * for heating in very high T models in map.in. zone 1 and 2 were 20% different, * since tau in is 1e-20, e2 is 0.9999, and so some H ots * these were not here at all*/ /* attenuation of flux by optical depths IN THIS ZONE * AngleIllum is 1/COS(theta), is usually 1, reset with illuminate command, * option for illumination of slab at an angle */ facexp.ExpZone[i] = sexp(opac.opac[i]*radius.drad/2.*Angllum.AngleIllum); /* e(-tau) in inward direction, up to illuminated face */ exptau.ExpmTau[i] = (float)facexp.ExpZone[i]; /* e2(tau) in inward direction, up to illuminated face */ e2tau.e2TauAbs[i] = (float)e2(deptau.depabs[i],exptau.ExpmTau[i]); rfield.reflin[i] = 0.; rfield.refcon[i] = 0.; rfield.reflux[i] = 0.; } /* update continuum */ SumContinuum(&SumOTS , &RatioOTSInc , &ipOTSchange); hydro.FreeFreeHeat = 0.; hydro.FreeFreeTotHeat = 0.; if( called.lgTalk ) { fprintf( ioQQQ, "1\n Start Iteration Number%2ld %75.75s\n\n\n", IterCnt.iter, input.chTitle ); } /* reset variable to do with FeII atom, in pop371 */ FeIIReset(); # ifdef DEBUG_FUN fputs( " <->RestartIter()\n", debug_fp ); # endif return; } /* do some work with ending the iteration */ void EndIter(void) { long i; float fac; # ifdef DEBUG_FUN fputs( "<+>RestartIter()\n", debug_fp ); # endif /* these were attenuated by too much in last zone */ for( i=0; i < NCELL; i++ ) { fac = (float)sexp(opac.opac[i]*radius.dReff/2.*Angllum.AngleIllum); if( fac > SMALLFLOAT ) { rfield.outcon[i] /= fac; rfield.outlin[i] /= fac; } } # ifdef DEBUG_FUN fputs( " <->EndIter()\n", debug_fp ); # endif return; }