/*EndZone last routine called after all zone calculations, before lgEndFun, * upon exit radiation field is for outer edge of current zone */ /*StartZone set variables that change with each zone, like radius, depth, * upon exit flux will be flux at center of zone about to be computed */ #include #include "cddefines.h" #include "phycon.h" #include "tlast.h" #include "opac.h" #include "rfield.h" #include "sign.h" #include "struc.h" #include "angllum.h" #include "sphere.h" #include "radius.h" #include "filfac.h" #include "globul.h" #include "tfidle.h" #include "pressure.h" #include "zonecnt.h" #include "tpred.h" #include "insane.h" #include "startendzone.h" /* this is number of zones to include in guess for next temperatures */ #define IOFF 3 void StartZone(char *chMode) { int lgNeOscil, lgTeOscil; long int i; double dt1 , de1, /* this is correction factor for dilution of radiation * across current zone, set in StartZone to correct * flux for center of current zone, and undone in ZoneDone */ DilutionCorrec , drFac , dTauThisZone, outwrd, ratio, rm, vin, vout; static double DTaver , DEaver, dt2 , de2; # ifdef DEBUG_FUN fputs( "<+>StartZone()\n", debug_fp ); # endif /* called at start of a zone to update all variables * having to do with starting this zone */ /* this sub can be called in two ways, 'incr' means increment * radius variables, while 'init' means only initialize rest */ if( strcmp(chMode,"init") == 0 ) { /* initialize the variables at start of calculation, after this exits * radius is equal to the initial radius, not outer edge of zone */ /* depth to illuminated face */ radius.depth = 1e-30; /* reset radius to inner radius, at this point equal to illuminated face radius */ radius.Radius = radius.rinner; /* thickness of next zone */ radius.drNext = radius.drad; /* depth variable for globule */ globul.glbdst = globul.glbrad; } else if( strcmp(chMode,"incr") == 0 ) { /* update radius variables */ radius.drad = radius.drNext; radius.depth += radius.drad; /* increment the radius, during the calculation Radius is the * outer radius of the current zone*/ radius.Radius += radius.drad*radius.dRadSign; /*********************************************************************** * * attenuuate rfield to center of this zone * ***********************************************************************/ /* radius was just incremented above, so this resets continuum to * flux at center of zone we are about to compute. radius will be * incremented immediately following this. this correction is undone * when ZoneDone called */ /* this will be the optical thickness of the next zone * AngleIllum is 1/COS(theta), is usually 1, reset with illuminate command, * option for illumination of slab at an angle */ /* radius.dRNeff is next dr with filling factor, this will only be * used to get approximate correction for attenuation * of radiation in this zone, * equations drived in hazy, netzer&ferland 83, for factor accounting * any changes here should be checked with both sphonly.in and pphonly*/ /* honlyotspp seems most sensitive to the 1.35 */ drFac = radius.drad*filfac.FillFac*Angllum.AngleIllum*1.35; /* dilute radiation so that rfield is in center of zone, this * goes into tmn at start of zone here, but is later divided out * when EndZone called */ DilutionCorrec = 1./POW2( (radius.Radius-radius.drad/2.)/(radius.Radius-radius.drad) ); /* note that this for loop is through <= nflux, to include the [nflux] * unit integration verification token. The token was set in * in MadeDiffuse and propagated out in metdif. the opacity at that energy is * zero, so only the geometrical part of tmn will affect things. The final * sum is verified in PrtComment */ for( i=0; i <= rfield.nflux; i++ ) { dTauThisZone = opac.opac[i]*drFac; /* TMN is array of scale factors which account for attenuation * of continuum across the zone (necessary to conserve energy * at the 1% - 5% level.) sphere effects in * drNext was set by NEXTDR and will be next dr */ if( dTauThisZone < 1e-4 ) { /* this small tau limit is the most likely branch, about 60% in parispn.in */ opac.tmn[i] = 1.f; } else if( dTauThisZone < 5. ) { /* this middle tau limit happens in the remaining 40% */ opac.tmn[i] = (float)((1. - exp(-dTauThisZone))/(dTauThisZone)); } else { /* this happens almost not at all */ opac.tmn[i] = (float)(1./(dTauThisZone)); } /* now add on correction for dilution across this zone */ opac.tmn[i] *= (float)DilutionCorrec; /* actually do the corrections */ rfield.flux[i] *= opac.tmn[i]; /*rfield.ConOutNoInter[i] *= opac.tmn[i]; rfield.outcon[i] *= opac.tmn[i]; rfield.outlin[i] *= opac.tmn[i];*/ } /* following is distance to globule, usually doesnt matter */ globul.glbdst -= (float)radius.drad; /* if gt 5th zone, and not constant pressure, and not oscillating te * then guess next temperature : option to predict next temperature * NZLIM is dim of structure variables saving temp, do data if nzone>NZLIM * IOFF set set to 5 in the define above */ if( (strcmp(pressure.chCPres,"CPRE") != 0 && ZoneCnt.nzone > IOFF + 1) && ZoneCnt.nzone < NZLIM ) { tpred.TeInit = phycon.te; tpred.EdenInit = (float)phycon.eden; lgTeOscil = FALSE; lgNeOscil = FALSE; dt1 = 0.; dt2 = 0.; de1 = 0.; de2 = 0.; DTaver = 0.; DEaver = 0.; for( i=ZoneCnt.nzone - IOFF-1; i < (ZoneCnt.nzone - 1); i++ ) { dt1 = dt2; de1 = de2; /* this will get the average change in temperature for the * past IOFF zones */ dt2 = struc.testr[i-1] - struc.testr[i]; de2 = struc.ednstr[i-1] - struc.ednstr[i]; DTaver += dt2; DEaver += de2; if( dt1*dt2 < 0. ) { lgTeOscil = TRUE; } if( de1*de2 < 0. ) { lgNeOscil = TRUE; } } /* option to guess next electron temperature if no oscillations */ if( !lgTeOscil ) { DTaver /= (float)(IOFF); /* don't want to overcorrect, use smaller of two */ dt2 = fabs(dt2); rm = fabs(DTaver); DTaver = sign(MIN2(rm,dt2),DTaver); /* do not let it change by more than 5% of current temp */ DTaver = sign(MIN2(rm,0.05*phycon.te),DTaver); /* this actually changes the temperature */ phycon.te -= (float)DTaver; } else { /* temp was oscillating - too dangerous to do anything */ DTaver = 0.; } /* this is used to remember the proposed temperature, for output * in the punch predictor command */ tpred.TeProp = phycon.te; /* option to guess next electron density if no oscillations */ if( !lgNeOscil ) { DEaver /= IOFF; de2 = fabs(de2); rm = fabs(DEaver); /* don't want to overcorrect, use smaller of two */ DEaver = sign(MIN2(rm,de2),DEaver); /* do not let it change by more than 5% of current temp */ DEaver = sign(MIN2(rm,0.05*phycon.eden),DEaver); /* this actually changes the temperature */ phycon.eden -= (float)DEaver; } else { /* temp was oscillating - too dangerous to do anything */ DEaver = 0.; } /* this is used to remember the proposed temperature, for output * in the punch predictor command */ tpred.EdenProp = (float)phycon.eden; tfidle(); } } else { fprintf( ioQQQ, " ZONSRT called with insane argument, %4.4s\n", chMode ); insane(); puts( "[Stop in StartZone]" ); exit(1); } /* check whether zone thickness is too small relative to radius */ if( strcmp(pressure.chCPres,"GLOB") == 0 ) { ratio = radius.drad/(globul.glbdst + radius.drad); /* this only matters for globule command */ if( ratio < 1e-14 ) { radius.lgdR2Small = TRUE; } else { radius.lgdR2Small = FALSE; } } filfac.FillFac = (float)(filfac.fiscal* pow( radius.Radius/radius.rinner ,filfac.filpow)); filfac.FillFac = (float)MIN2(1.,filfac.FillFac); /* effective radius */ radius.dReff = radius.drad*filfac.FillFac; /* area factor, ratio of radius to out edge of this zone * relative to radius of illuminated face of cloud */ /*radius.r1r0sq = (float)POW2( (radius.Radius - radius.drad*radius.dRadSign/2.)/radius.rinner);*/ /*>>>chng 99 jan 25, definition of r1r0sq is now outer edge of current zone * relative to inner edge of cloud */ radius.r1r0sq = POW2( radius.Radius/radius.rinner ); /* following only approximate, used for center of next zone */ radius.dRNeff = radius.drNext*filfac.FillFac; /* this is the middle of the zone */ rm = radius.Radius - radius.drad/2.; /* RADIUS is outer radius of this zone, so radius - drad is inner radius * RINNER is inner radius of nebula; dVeff is the effective vol element * dVeff is vol rel to inner radius, so has units cm * for plane parallel dVeff is dReff */ if( radius.drad/radius.Radius > 0.01 ) { vin = ((radius.Radius - radius.drad)/radius.rinner)* (MIN2((radius.Radius-radius.drad),radius.CylindHigh)/radius.rinner)* (radius.Radius - radius.drad); vout = (radius.Radius/radius.rinner)* (MIN2(radius.Radius,radius.CylindHigh)/radius.rinner)* radius.Radius; /* this is possible correction for slit projected onto resolved object - * this only happens when sphere slit is entered. * NB not clear what happens when slit and cylendar are both used - * this would be very orientation specific */ /* just the difference in the two vols */ /*radius.dVeff = (vout - vin)/3.*filfac.FillFac;*/ radius.dVeff = (vout - vin)/3.*filfac.FillFac; if( sphere.lgSlit ) { /* default of iEmissPower is 2, set to 0 is just sphere slit, * and to 1 if sphere long slit set */ if( sphere.iEmissPower == 2 ) { /* this is long slit option, remove one power of radius */ radius.dVeff /= (rm/radius.rinner); } else if( sphere.iEmissPower == 1 ) { /* this is short slit option, remove two powers of radius */ radius.dVeff /= POW2(rm/radius.rinner); } } } else { /* RM is the middle of the zone; * radius is always the OUTER edge of the current zone */ radius.dVeff = (rm/radius.rinner)*radius.drad*filfac.FillFac* (MIN2(rm,radius.CylindHigh)/radius.rinner); /* slit option */ if( sphere.lgSlit ) { /* default of iEmissPower is 3, set to 1 is just sphere beam, * and to 2 if sphere long slit set */ if( sphere.iEmissPower == 2 ) { /* this is long slit option, remove one power of radius */ radius.dVeff /= (rm/radius.rinner); } else if( sphere.iEmissPower == 1 ) { /* this is short slit option, remove two powers of radius */ radius.dVeff /= POW2(rm/radius.rinner); } } } /* covering factor, default is unity */ radius.dVeff *= sphere.covgeo; /* these are needed for line intensities in outward beam * lgSphere set, sphere.covrt usually 1, 0 when not lgSphere * so outwrd is 1 when lgSphere set 1/2 when not set, */ outwrd = (1. + sphere.covrt)/2.; /*>>>chng 99 apr 23, from above to below */ radius.dVolOutwrd = outwrd*POW2( (radius.Radius-radius.dReff/2.)/radius.Radius) * radius.drad; /* this includes covering fact, the effective vol,, and 1/r^2 dilution, * dReff includes filling factor. the covering factor part is 1 for sphere, * 1/2 for open */ radius.dVolOutwrd = outwrd*radius.Radius*radius.dReff/(radius.Radius + 2.*radius.drad); /* dReff from above, includes filling factor */ radius.dVolOutwrd = outwrd*radius.dReff; /* following will be used to "uncorrect" for this in lines when predicting continua*/ radius.GeoDil = radius.Radius/(radius.Radius + 2.*radius.drad); /* this should multiply the line to become the net inward. geo part is 1/2 for * open geomery, 0 for closed. for case of isotropic emitter only this and dVolOutwrd * above are used */ radius.dVolReflec = (1. - outwrd)*radius.dReff*radius.r1r0sq; if( sphere.lgSphere ) { /* inward beams do not go in when lgSphere set since geometry symmetric */ radius.BeamInIn = 0.; radius.BeamInOut = radius.Radius*radius.dReff/(radius.Radius + 2.*radius.drad); } else { radius.BeamInIn = radius.dReff*radius.r1r0sq; /* inward beams do not go out */ radius.BeamInOut = 0.; } /* factor for outwardly directed part of line */ radius.BeamOutOut = radius.Radius*radius.dReff/(radius.Radius + 2.*radius.drad); # ifdef DEBUG_FUN fputs( " <->StartZone()\n", debug_fp ); # endif return; } void EndZone(void) { long i; # ifdef DEBUG_FUN fputs( "<+>EndZone()\n", debug_fp ); # endif /* this is the last routine called after all zone calculations, * and before calling lgEndFun to check if done */ tlastCom.tlast = phycon.te; /*********************************************************************** * * correct rfield for attenuation from center of zone to inner edge * ***********************************************************************/ /* radius is outer radius of this zone, this resets continuum to * flux at illuminated face of zone we have already computed. */ /* opac.tmn defined in StartZone */ /* undilute radiation so that rfield is at face of zone */ /* NB, upper limit of sum includes [nflux] to test unit verification cell */ for( i=0; i <= rfield.nflux; i++ ) { rfield.flux[i] /= opac.tmn[i]; /*rfield.ConOutNoInter[i] /= opac.tmn[i]; rfield.outcon[i] /= opac.tmn[i]; rfield.outlin[i] /= opac.tmn[i];*/ } # ifdef DEBUG_FUN fputs( " <->EndZone()\n", debug_fp ); # endif return; }