/*ConvInitTemp drive search for initial temperature, for illuminated face */ #include #include #include "cddefines.h" #include "trace.h" #include "mappar.h" #include "elecden.h" #include "stopcalc.h" #include "struc.h" #include "nodec.h" #include "addopac.h" #include "colneg.h" #include "rfield.h" #include "thigh.h" #include "heavy.h" #include "forcte.h" #include "h.h" #include "densty.h" #include "radius.h" #include "tcool.h" #include "itercnt.h" #include "phycon.h" #include "tfidle.h" #include "boltgn.h" #include "heat.h" #include "pressure.h" #include "showme.h" #include "hevmolec.h" #include "ionfracs.h" #include "punt.h" #include "converge.h" /* XMULTP is scale factor to multiply temperature */ #define XMULTP 3.0 int ConvInitTemp(void) { long int i, ionlup, isear, j, looplimit, nelem , nflux_old, nelem_reflux , ion_reflux; double cnew, cold, dt, fac2, factor, hnew, hold, sign, t1, t2, tinc, tmax, tmin, tnew, told; static double TempPrevIteration = -DBL_MAX; double derrdt; # ifdef DEBUG_FUN fputs( "<+>ConvInitTemp()\n", debug_fp ); # endif tnew = DBL_MAX; cnew = DBL_MAX; hnew = DBL_MAX; /* this counts number of times ionize is called by PresChng, in current zone */ conv.nPres2Ioniz = 0; /* this counts how many times ionize is called in this model after startr, * and is flag used by ionize to understand it is being called the first time*/ conv.nTotalIoniz = 0; /* ots rates not oscillating */ conv.lgOscilOTS = FALSE; ElecDen.lgEdenBad = FALSE; ElecDen.nzEdenBad = 0; /* these are variables to remember the biggest error in the * electron density, and the zone in qhich it occurred */ conv.BigEdenError = 0.; conv.AverEdenError = 0.; conv.BigHeatCoolError = 0.; conv.AverHeatCoolError = 0.; conv.BigPressError = 0.; conv.AverPressError = 0.; if( trace.lgTrConvg ) { fprintf( ioQQQ, " \n" ); fprintf( ioQQQ, " ConvInitTemp entered \n" ); } # if 0 /* as a test never do map from here, since caught later, * must also put else in following if if turned back on */ /******************************************************************** * * if MapZone LT 0, just map first zone and punch out * *********************************************************************/ if( MapPar.MapZone <= 0 ) { FILE *ioSAVE; if( trace.lgTrace || trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp called. just map\n" ); } ioSAVE = ioQQQ; ioQQQ = ioMAP; punt(" map"); ioQQQ = ioSAVE; } # endif /******************************************************************** * * force an initial temperature, also constant temperature * *********************************************************************/ if( forcte.ForceTemp != 0. ) { /* force te to certain value */ if( trace.lgTrace ) { fprintf( ioQQQ, " ConvInitTemp called, forcing TE to%10.2e\n", forcte.ForceTemp ); } phycon.te = forcte.ForceTemp; /* evaluate current pressure for this constant temperature model */ pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); /* say this is search for first solution, let Ne change a lot */ conv.lgSearch = TRUE; ionlup = 0; conv.lgConvIoniz = FALSE; /* >>chng 96 oct 11, loop at least 3x to make sure ionization converged */ while( ionlup <= 3 || (!conv.lgConvIoniz && ionlup <= 7) ) { /* update density - temperature variables which may have changed */ tfidle(); /* reevaluate gaunt factors and Boltzmann factors */ boltgn(); /* lgIonDone can be set false by several tests in ionize */ conv.lgConvIoniz = TRUE; ConvEdenIoniz(); /* only used to keep forcing new optical depths below */ ionlup += 1; } if( trace.lgTrace || trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp: search set false\n" ); } conv.lgSearch = FALSE; } /******************************************************************** * * this is second or higher iteration, reestablish original temperature * *********************************************************************/ else if( IterCnt.iter != 1 ) { /* this is second or higher iteration on multi-iteration model */ if( trace.lgTrace || trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp called, ITER=%2ld resetting Te to %10.2e\n", IterCnt.iter, TempPrevIteration ); } if( trace.lgTrace || trace.lgTrConvg ) { fprintf( ioQQQ, " search set true\n" ); } /* search phase must be turned on so that variables such as the ors rates, * secondary ionization, and auger yields, can converge more quickly to * proper values */ conv.lgSearch = TRUE; /* this is the temperature and pressure from the previous iteration */ phycon.te = (float)TempPrevIteration; /* the inital pressure should now be valid */ pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); /* this counts how many times ionize is called in this model after startr, * and is flag used by ionize to understand it is being called the first time*/ conv.nTotalIoniz = 0; /* now get final temperature */ ConvTempEdenIoniz(); if( trace.lgTrace || trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp: search set false\n" ); } conv.lgSearch = FALSE; /* now get final temperature */ ConvTempEdenIoniz(); /* the inital pressure should now be valid */ pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); /* remember the current temperature for next time */ TempPrevIteration = phycon.te; } else { /******************************************************************** * * do first te from scratch * *********************************************************************/ /* say that we are in search phase */ conv.lgSearch = TRUE; if( trace.lgTrace ) { fprintf( ioQQQ, " ConvInitTemp called, new temperature.\n" ); } /* approach equilibrium from either high or low TE */ if( thigh.lgTeHigh ) { /* approach from high TE */ phycon.te = 1e6; nodec.lgNoDec = FALSE; factor = 1./XMULTP; } else { /* approach from low TE */ nodec.lgNoDec = TRUE; phycon.te = 4000.; factor = XMULTP; } /* find first temperature */ colneg.lgColNeg = FALSE; pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); heat.htot = 1.; tcool.ctot = 1.; /* call cooling, heating, opacity, loop to convergence * this is very first call to it, by default is at 4000K */ if( trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp calling ConvEdenIoniz1 with te=%10.3e\n", phycon.te ); } /* update density - temperature variables which may have changed */ tfidle(); /* reevaluate gaunt factors and Boltzmann factors */ boltgn(); ConvEdenIoniz(); /* >>chng 97 feb 18, define following three */ told = phycon.te; cold = tcool.ctot; hold = heat.htot; if( trace.lgTrace ) { fprintf( ioQQQ, " RETURN TO ConvInitTemp 1, HTOT:%10.3e CTOT:%10.3e Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", heat.htot, tcool.ctot, phycon.te, phycon.eden ); } /* sign is cooling minus heating, in most cases negative means T too low */ sign = fabs(tcool.ctot-heat.htot)/(tcool.ctot - heat.htot); /* if low te approach then check if already past soln (cooling > heating) */ if( !thigh.lgTeHigh && sign >= 0. ) { if( trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp entering T>chng 96 may 30, logic had gone straight to 3K, which blew up dense clouds * now go down more gently * >>chng 97 mar 3, variable factor, smaller when close to 3K */ fac2 = 0.8; while( phycon.te*fac2 > StopCalc.TeLowest && sign >= 0. ) { /* >>chng 96 June 14, to below, so that we can get soln near 3K * do while( te.gt.TeLowest .and. sign.ge.0. ) * >>chng 97 feb 18, define following three */ told = phycon.te; cold = tcool.ctot; hold = heat.htot; /* >>chng 97 mar 3, make fac2 small when close to lower bound */ if( phycon.te < 5. ) { fac2 = 0.97; } else if( phycon.te < 20. ) { fac2 = 0.9; } else { fac2 = 0.8; } /* now change to new temp */ phycon.te *= (float)fac2; /* te = TeLowest */ pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); /* call cooling, heating, opacity, loop to convergence */ if( trace.lgTrConvg ) { fprintf( ioQQQ, "\n calling ConvEdenIoniz2, te=%9.2e\n", phycon.te ); } /* update density - temperature variables which may have changed */ tfidle(); /* reevaluate gaunt factors and Boltzmann factors */ boltgn(); ConvEdenIoniz(); tnew = phycon.te; cnew = tcool.ctot; hnew = heat.htot; if( trace.lgTrConvg ) { fprintf( ioQQQ, " ConvEdenIoniz2 returns, te=%9.2e htot%9.2e ctot%9.2e\n", phycon.te, heat.htot, tcool.ctot ); } if( trace.lgTrace || trace.lgTrConvg ) { fprintf( ioQQQ, " RETURN TO ConvInitTemp 2 HTOT:%10.3e CTOT:%10.3e Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", heat.htot, tcool.ctot, phycon.te, phycon.eden ); } sign = fabs(tcool.ctot-heat.htot)/(tcool.ctot - heat.htot); } if( tcool.ctot > heat.htot ) { /* even at lowest temp, too much cooling; no solution possible */ fprintf( ioQQQ, " Equilibrium temperture below T=%8.2e. If this is OK, then reset with LOWEST command. This is ConvInitTemp. TeLowest was%9.2e\n", phycon.te, StopCalc.TeLowest ); # ifdef DEBUG_FUN fputs( " <->ConvInitTemp()\n", debug_fp ); # endif /* this is an error return, calculation will immediately stop */ return(1); } /* >>chng 97 mar 3, from *0.9 to / fac2 */ phycon.te /= (float)fac2; factor = 1./(0.5*(1. + fac2)); sign = fabs(cold-hold)/(cold - hold); if( trace.lgTrConvg ) { fprintf( ioQQQ, " <4000K loop returns, te=%9.2e sign%9.2e factor%9.2e\n", phycon.te, sign, factor ); } } tmin = 1e9; tmax = 0.; /* coming up to here te is either 4,000 (usually) or 10^6 */ t1 = StopCalc.TeLowest/1.001; /* >>chng 96 dc 28, added div by factor to prevent checking on temp > 1e10 */ t2 = StopCalc.TeHighest*1.001/factor; if( trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp entering second loop te=%9.2e sign%9.2e factor%9.2e t1, t2=%10.2e%10.2e\n", phycon.te, sign, factor, t1, t2 ); } while( ((tcool.ctot - heat.htot)*sign > 0. && phycon.te > t1) && phycon.te < t2 ) { /* >>chng 97 feb 18, define following three */ told = phycon.te; cold = tcool.ctot; hold = heat.htot; phycon.te *= (float)factor; pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); tmin = MIN2(tmin,phycon.te); tmax = MAX2(tmax,phycon.te); conv.lgConvIoniz = FALSE; strcpy( conv.chConvIoniz, "ConvInitTemp srh" ); isear = 1; while( isear <= 6 && !conv.lgConvIoniz ) { /* call cooling, heating, opacity, loop to convergence */ conv.lgConvIoniz = TRUE; if( trace.lgTrConvg ) { fprintf( ioQQQ, "\n calling ConvEdenIoniz3, te=%9.2e\n", phycon.te ); } /* update density - temperature variables which may have changed */ tfidle(); /* reevaluate gaunt factors and Boltzmann factors */ boltgn(); ConvEdenIoniz(); if( trace.lgTrConvg ) { fprintf( ioQQQ, " ConvEdenIoniz3 returns, te=%9.2e htot%9.2e ctot%9.2e\n", phycon.te, heat.htot, tcool.ctot ); } if( trace.lgTrace ) { fprintf( ioQQQ, " RETURN TO ConvInitTemp 3 HTOT:%10.3e ctot:%10.3e Te:%11.4e eden%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", heat.htot, tcool.ctot, phycon.te, phycon.eden ); } isear += 1; } tnew = phycon.te; cnew = tcool.ctot; hnew = heat.htot; } /* fall through loop, te too high or low */ if( phycon.te <= t1 || phycon.te >= t2 ) { if( phycon.te <= t1 ) { fprintf( ioQQQ, " ConvInitTemp finds TE below TE = %10.3eK, HTOT=%12.4e CTOT=%12.4e\n", t1, heat.htot, tcool.ctot ); } else if( phycon.te >= t2 ) { fprintf( ioQQQ, " ConvInitTemp finds TE above TE = %10.3eK, HTOT=%12.4e last CTOT=%12.4e\n", t2, heat.htot, tcool.ctot ); } else { fprintf( ioQQQ, " ConvInitTemp finds insanity \n" ); ShowMe(); } # ifdef DEBUG_FUN fputs( " <->ConvInitTemp()\n", debug_fp ); # endif /* this is an error condition, calculation will immediately stop */ return(1); } /* now search inside (also either side of) TE range where sign changed * first turn on ionization stage trimming */ nodec.lgNoDec = FALSE; colneg.lgColNeg = TRUE; /* change sense of temp search so that we can "weave" into soln */ tinc = 1./pow(factor,0.2); for( j=0; j < 2; j++ ) { if( (tcool.ctot - heat.htot > 1e-30) && tcool.ctot > 0. ) { sign = fabs(tcool.ctot-heat.htot)/(tcool.ctot - heat.htot); if( trace.lgTrConvg ) { fprintf( ioQQQ, "\n\n >>>entering second loop, te reset to%9.2e\n", phycon.te ); } pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); /* call cooling, heating, opacity, loop to convergence */ isear = 1; conv.lgConvIoniz = FALSE; if( trace.lgTrace ) { fprintf( ioQQQ, " RETURN TO ConvInitTemp 4 HTOT:%10.3e CTOT:%10.3e Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", heat.htot, tcool.ctot, phycon.te, phycon.eden ); } i = 0; /* >>chng 97 july 29, to recover fumbled initial ionization */ looplimit = 12; told = phycon.te; cold = tcool.ctot; hold = heat.htot; /* this loop is not necessariy executed (ctot and htot very close) */ while( i < looplimit && (tcool.ctot - heat.htot)*sign > 0. ) { i += 1; told = phycon.te; cold = tcool.ctot; hold = heat.htot; phycon.te *= (float)tinc; pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); /* update density - temperature variables which may have changed */ tfidle(); /* reevaluate gaunt factors and Boltzmann factors */ boltgn(); ConvEdenIoniz(); if( trace.lgTrConvg ) { fprintf( ioQQQ, " ConvEdenIoniz returned, te=%9.2e htot%9.2e ctot%9.2e\n", phycon.te, heat.htot, tcool.ctot ); } tnew = phycon.te; cnew = tcool.ctot; hnew = heat.htot; if( trace.lgTrace ) { fprintf( ioQQQ, " RETURN TO ConvInitTemp 4 HTOT:%10.3e CTOT:%10.3e Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", heat.htot, tcool.ctot, phycon.te, phycon.eden ); } } if( (i >= looplimit && (tcool.ctot - heat.htot)*sign > 0.) && fabs(tcool.ctot/heat.htot-1.) > 0.05 ) { /* >>chng 96 dec 17, way out for dense, cold, molecular gas */ if( phycon.te < 1e3 && hevmolec.hevmol[IPCO-1]/xIonFracs[0][5] > 0.5 ) { fprintf( ioQQQ, " cold molecular trap, te=%10.2e F(CO)=%10.2e\n", phycon.te, hevmolec.hevmol[IPCO-1]/xIonFracs[0][5] ); fprintf( ioQQQ, " Conditions in this cloud are not appropriate for this code.\n" ); # ifdef DEBUG_FUN fputs( " <->ConvInitTemp()\n", debug_fp ); # endif /* this is an error return, calculation will immediately stop */ return(1); } else { fprintf( ioQQQ, " Insanity occurred in ConvInitTemp, i=0\n" ); ShowMe(); puts( "[Stop in ConvInitTemp]" ); exit(1); } } /* weave back in again */ tinc = 1./pow(tinc,0.2); } } if( trace.lgTrace || trace.lgTrConvg ) { fprintf( ioQQQ, " ConvInitTemp: search set false\n" ); } conv.lgSearch = FALSE; /* now must be close enough to call regular ionization-cooling * first do linear interpolation to guess temperature */ if( (tnew != told && tnew > 0.) && told > 0. ) { derrdt = (double)((cnew-hnew)-(cold-hold))/(double)(tnew-told); dt = -(double)(cold-hold)/derrdt; phycon.te = (float)(told + dt); phycon.te = (float)MAX2(phycon.te,MIN2(told,tnew)); phycon.te = (float)MIN2(phycon.te,MAX2(told,tnew)); pressure.PressureInit = (float)TotalPressure(); pressure.RamPresInit = (float)RamPressure(); } if( trace.lgTrace ) { fprintf( ioQQQ, " ConvInitTemp calls IONTE after interpolating a temperature of%10.3e Told=%10.2e Tnew=%10.2e dC-Hold%10.2e dC-Hnew%10.2e\n", phycon.te, told, tnew, cold - hold, cnew - hnew ); } if( trace.lgTrConvg ) { fprintf( ioQQQ, "\n\n ConvInitTemp calling ConvTempEdenIoniz with interpolated te=%10.3e\n", phycon.te ); } ConvTempEdenIoniz(); TempPrevIteration = phycon.te; if( trace.lgTrace ) { fprintf( ioQQQ, " ConvInitTemp return, TE:%10.2e==================\n", phycon.te ); } } /* this counts number of times ionize is called by PresChng, in current zone * these are reset here, so that we count from first zone not search */ conv.nPres2Ioniz = 0; /* this counts how many times ionize is called in this model after startr, * and is flag used by ionize to understand it is being called the first time*/ conv.nTotalIoniz = 0; ElecDen.lgEdenBad = FALSE; ElecDen.nzEdenBad = 0; /* these are variables to remember the biggest error in the * electron density, and the zone in qhich it occurred */ conv.BigEdenError = 0.; conv.AverEdenError = 0.; conv.BigHeatCoolError = 0.; conv.AverHeatCoolError = 0.; conv.BigPressError = 0.; conv.AverPressError = 0.; /* now remember some things we may need even in first zone, these * are normally set towards end of zone calc in tauinc */ struc.testr[0] = phycon.te; struc.pdenstr[0] = densty.pden; struc.heatstr[0] = heat.htot; struc.coolstr[0] = tcool.ctot; struc.volstr[0] = (float)radius.dVeff; struc.radstr[0] = (float)radius.dReff; struc.histr[0] = h.hi; struc.hiistr[0] = h.hii; struc.ednstr[0] = (float)phycon.eden; struc.o3str[0] = xIonFracs[3][7]; struc.DenMass[0] = densty.xMassDensity; /* check that nflux extends above IP of highest ionization species present. * for collisional case it is possible for species to exist that are higher * IP than the limit to the continuum. Need continuum to encompass all * possible emission */ nflux_old = rfield.nflux; nelem_reflux = -1; ion_reflux = -1; for( nelem=2; nelem < LIMELM; nelem++ ) { /* do not include hydrogenic species in following */ for( i=0; i < nelem; i++ ) { if( xIonFracs[i+2][nelem] > 0. ) { if( Heavy.ipHeavy[i][nelem] > rfield.nflux ) { rfield.nflux = Heavy.ipHeavy[i][nelem] ; nelem_reflux = nelem; ion_reflux = i; } } } } /* was the upper limit to the opacities updated? */ if( nflux_old != rfield.nflux ) { /* update opacities, can't have zeros across band between nflux_old and nflux */ AddOpac(); /* need to update continuum opacities */ if( trace.lgTrace ) { fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n", nflux_old , rfield.nflux , rfield.anu[nflux_old] , rfield.anu[rfield.nflux] ); fprintf(ioQQQ," caused by element %li ion %li \n", nelem_reflux ,ion_reflux ); } } # ifdef DEBUG_FUN fputs( " <->ConvInitTemp()\n", debug_fp ); # endif /* this is the only valid return */ return(0); }