00001 /* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland 00002 * For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ZoneEnd last routine called after all zone calculations, before iter_end_check, 00004 * upon exit radiation field is for outer edge of current zone */ 00005 /*ZoneStart set variables that change with each zone, like radius, depth, 00006 * upon exit flux will be flux at center of zone about to be computed */ 00007 #include "cddefines.h" 00008 #include "phycon.h" 00009 #include "opacity.h" 00010 #include "rfield.h" 00011 #include "struc.h" 00012 #include "thermal.h" 00013 #include "dense.h" 00014 #include "h2.h" 00015 #include "geometry.h" 00016 #include "conv.h" 00017 #include "dynamics.h" 00018 #include "radius.h" 00019 #include "zones.h" 00020 /* this is number of zones to include in guess for next temperatures */ 00021 #define IOFF 3 00022 00023 void ZoneStart(const char *chMode) 00024 { 00025 bool lgNeOscil, 00026 lgTeOscil; 00027 long int i; 00028 00029 double dt1 , de1, 00030 /* this is correction factor for dilution of radiation 00031 * across current zone, set in ZoneStart to correct 00032 * flux for center of current zone, and undone in ZoneDone */ 00033 DilutionCorrec , 00034 drFac , 00035 dTauThisZone, 00036 outwrd, 00037 ratio, 00038 rm, 00039 rad_middle_zone, 00040 vin, 00041 vout; 00042 00043 static double DTaver , DEaver, 00044 dt2 , de2; 00045 00046 DEBUG_ENTRY( "ZoneStart()" ); 00047 00048 /* this sub can be called in two ways, 'incr' means increment 00049 * radius variables, while 'init' means only initialize rest */ 00050 /* called at start of a zone to update all variables 00051 * having to do with starting this zone */ 00052 00053 /* first establish current filling factor (normally 1) since used within 00054 * following branches */ 00055 geometry.FillFac = (float)(geometry.fiscal* 00056 pow( radius.Radius/radius.rinner ,(double)geometry.filpow)); 00057 geometry.FillFac = (float)MIN2(1.,geometry.FillFac); 00058 00059 if( strcmp(chMode,"init") == 0 ) 00060 { 00061 /* initialize the variables at start of calculation, after this exits 00062 * radius is equal to the initial radius, not outer edge of zone */ 00063 00064 /* depth to illuminated face */ 00065 radius.depth = 1e-30; 00066 00067 /* integral of depth times filling factor */ 00068 radius.depth_x_fillfac = radius.depth*geometry.FillFac; 00069 /* effective radius */ 00070 radius.drad_x_fillfac = radius.drad*geometry.FillFac; 00071 00072 /* reset radius to inner radius, at this point equal to illuminated face radius */ 00073 radius.Radius = radius.rinner; 00074 radius.Radius_mid_zone = radius.rinner + radius.drad/2.; 00075 00076 /* thickness of next zone */ 00077 radius.drNext = radius.drad; 00078 00079 00080 /* depth to illuminated face */ 00081 radius.depth_mid_zone = radius.drad/2.; 00082 00083 /* depth variable for globule */ 00084 radius.glbdst = radius.glbrad; 00085 00086 /* this is case where outer radius is set */ 00087 if( radius.router[iteration-1] < 5e29 ) 00088 { 00089 radius.Depth2Go = radius.router[iteration-1]; 00090 } 00091 else if( iteration > 1 ) 00092 { 00093 /* this case second or later iteration, use previous thickness */ 00094 radius.Depth2Go = radius.router[iteration-2]; 00095 } 00096 else 00097 { 00098 /* first iteration and depth not set, so we cannot estimate it */ 00099 radius.Depth2Go = 0.; 00100 } 00101 } 00102 else if( strcmp(chMode,"incr") == 0 ) 00103 { 00104 /* update radius variables - called by cloudy at start of this zone's calcs */ 00105 radius.drad = radius.drNext; 00106 /* >>chng 06 mar 21, remember largest and smallest dr over this iteration 00107 * for use in recombination time dep logic */ 00108 radius.dr_min_last_iter = MIN2( radius.dr_min_last_iter , radius.drNext ); 00109 radius.dr_max_last_iter = MAX2( radius.dr_max_last_iter , radius.drNext ); 00110 00111 ASSERT( radius.drad > 0. ); 00112 radius.depth += radius.drad; 00113 radius.depth_mid_zone = radius.depth - radius.drad/2.; 00114 00115 /* effective radius */ 00116 radius.drad_x_fillfac = radius.drad*geometry.FillFac; 00117 00118 /* integral of depth times filling factor */ 00119 radius.depth_x_fillfac += radius.drad_x_fillfac; 00120 00121 /* decrement Depth2Go but do not let become negative */ 00122 radius.Depth2Go = MAX2( 0.,radius.Depth2Go - radius.drad ); 00123 00124 /* increment the radius, during the calculation Radius is the 00125 * outer radius of the current zone*/ 00126 radius.Radius += radius.drad*radius.dRadSign; 00127 00128 /* Radius is now outer edge of this zone since incremented above, 00129 * so need to add drad to it */ 00130 radius.Radius_mid_zone = radius.Radius - radius.drad/2.; 00131 00132 /*********************************************************************** 00133 * 00134 * attenuate rfield to center of this zone 00135 * 00136 ***********************************************************************/ 00137 00138 /* radius was just incremented above, so this resets continuum to 00139 * flux at center of zone we are about to compute. radius will be 00140 * incremented immediately following this. this correction is undone 00141 * when ZoneDone called */ 00142 00143 /* this will be the optical thickness of the next zone 00144 * AngleIllumRadian is 1/COS(theta), is usually 1, reset with illuminate command, 00145 * option for illumination of slab at an angle */ 00146 00147 /* radius.dRNeff is next dr with filling factor, this will only be 00148 * used to get approximate correction for attenuation 00149 * of radiation in this zone, 00150 * equations derived in hazy, Netzer&Ferland 83, for factor accounting 00151 * any changes here should be checked with both sphonly.in and pphonly*/ 00152 /* honlyotspp seems most sensitive to the 1.35 */ 00153 drFac = radius.drad*geometry.FillFac*geometry.DirectionalCosin*1.35; 00154 00155 /* dilute radiation so that rfield is in center of zone, this 00156 * goes into tmn at start of zone here, but is later divided out 00157 * when ZoneEnd called */ 00158 DilutionCorrec = 1./POW2( 00159 (radius.Radius-radius.drad/2.)/(radius.Radius-radius.drad) ); 00160 00161 /* note that this for loop is through <= nflux, to include the [nflux] 00162 * unit integration verification token. The token was set in 00163 * in MadeDiffuse and propagated out in metdif. the opacity at that energy is 00164 * zero, so only the geometrical part of tmn will affect things. The final 00165 * sum is verified in PrtComment */ 00166 for( i=0; i <= rfield.nflux; i++ ) 00167 { 00168 dTauThisZone = opac.opacity_abs[i]*drFac; 00169 /* TMN is array of scale factors which account for attenuation 00170 * of continuum across the zone (necessary to conserve energy 00171 * at the 1% - 5% level.) sphere effects in 00172 * drNext was set by NEXTDR and will be next dr */ 00173 00174 if( dTauThisZone < 1e-4 ) 00175 { 00176 /* this small tau limit is the most likely branch, about 60% in parispn.in */ 00177 opac.tmn[i] = 1.f; 00178 } 00179 else if( dTauThisZone < 5. ) 00180 { 00181 /* this middle tau limit happens in the remaining 40% */ 00182 opac.tmn[i] = (float)((1. - exp(-dTauThisZone))/(dTauThisZone)); 00183 } 00184 else 00185 { 00186 /* this happens almost not at all */ 00187 opac.tmn[i] = (float)(1./(dTauThisZone)); 00188 } 00189 00190 /* now add on correction for dilution across this zone */ 00191 opac.tmn[i] *= (float)DilutionCorrec; 00192 00193 /* actually do the corrections */ 00194 rfield.flux[i] *= opac.tmn[i]; 00195 00196 /* >>chng 03 nov 08, update SummedCon here since flux changed */ 00197 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i]; 00198 } 00199 00200 /* following is distance to globule, usually does not matter */ 00201 radius.glbdst -= (float)radius.drad; 00202 00203 /* if gt 5th zone, and not constant pressure, and not oscillating te 00204 * then guess next temperature : option to predict next temperature 00205 * NZLIM is dim of structure variables saving temp, do data if nzone>NZLIM 00206 * IOFF is number of zones to look over, set set to 3 in the define above */ 00207 /* >>chng 01 mar 12, add option to not do this, set with no tepred command */ 00208 if( (strcmp(dense.chDenseLaw,"CPRE") != 0) && 00209 thermal.lgPredNextTe && 00210 (nzone > IOFF + 1) ) 00211 { 00212 phycon.TeInit = phycon.te; 00213 phycon.EdenInit = dense.eden; 00214 lgTeOscil = false; 00215 lgNeOscil = false; 00216 dt1 = 0.; 00217 dt2 = 0.; 00218 de1 = 0.; 00219 de2 = 0.; 00220 DTaver = 0.; 00221 DEaver = 0.; 00222 for( i=nzone - IOFF-1; i < (nzone - 1); i++ ) 00223 { 00224 dt1 = dt2; 00225 de1 = de2; 00226 /* this will get the average change in temperature for the 00227 * past IOFF zones */ 00228 dt2 = struc.testr[i-1] - struc.testr[i]; 00229 de2 = struc.ednstr[i-1] - struc.ednstr[i]; 00230 DTaver += dt2; 00231 DEaver += de2; 00232 if( dt1*dt2 < 0. ) 00233 { 00234 lgTeOscil = true; 00235 } 00236 if( de1*de2 < 0. ) 00237 { 00238 lgNeOscil = true; 00239 } 00240 } 00241 00242 /* option to guess next electron temperature if no oscillations */ 00243 if( !lgTeOscil ) 00244 { 00245 DTaver /= (double)(IOFF); 00246 /* don't want to over correct, use smaller of two */ 00247 dt2 = fabs(dt2); 00248 rm = fabs(DTaver); 00249 DTaver = sign(MIN2(rm,dt2),DTaver); 00250 /* do not let it change by more than 5% of current temp */ 00251 /* >>chng 03 mar 18, from 5% to 1% - convergence is much better 00252 * now, and throwing the next Te off by 5% could easily disturb 00253 * the solvers - primal.in was a good case */ 00254 DTaver = sign(MIN2(rm,0.01*phycon.te),DTaver); 00255 /* this actually changes the temperature */ 00256 phycon.te = phycon.te - DTaver; 00257 tfidle(true); 00258 } 00259 else 00260 { 00261 /* temp was oscillating - too dangerous to do anything */ 00262 DTaver = 0.; 00263 } 00264 /* this is used to remember the proposed temperature, for output 00265 * in the punch predictor command */ 00266 phycon.TeProp = phycon.te; 00267 00268 /* option to guess next electron density if no oscillations */ 00269 if( !lgNeOscil ) 00270 { 00271 DEaver /= IOFF; 00272 de2 = fabs(de2); 00273 rm = fabs(DEaver); 00274 /* don't want to over correct, use smaller of two */ 00275 DEaver = sign(MIN2(rm,de2),DEaver); 00276 /* do not let it change by more than 5% of current temp */ 00277 DEaver = sign(MIN2(rm,0.05*dense.eden),DEaver); 00278 /* this actually changes the temperature */ 00279 dense.eden -= DEaver; 00280 } 00281 else 00282 { 00283 /* temp was oscillating - too dangerous to do anything */ 00284 DEaver = 0.; 00285 } 00286 /* this is used to remember the proposed temperature, for output 00287 * in the punch predictor command */ 00288 phycon.EdenProp = dense.eden; 00289 tfidle(false); 00290 } 00291 } 00292 00293 else 00294 { 00295 fprintf( ioQQQ, " PROBLEM ZoneStart called with insane argument, %4.4s\n", 00296 chMode ); 00297 /* TotalInsanity exits after announcing a problem */ 00298 TotalInsanity(); 00299 } 00300 00301 /* find the mean dr for the past few zones - this var is used 00302 * in setting the optical scale in induced processes, and large 00303 * changes in it will drive oscillations in the H+ fraction */ 00304 radius.drad_x_fillfac_mean = radius.drad_x_fillfac; 00305 i = 1; 00306 while( i < 5 && nzone-i-1>=0 ) 00307 { 00308 radius.drad_x_fillfac_mean += struc.drad_x_fillfac[nzone-i-1]; 00309 ++i; 00310 } 00311 radius.drad_x_fillfac_mean /= i; 00312 00313 /* do advection if enabled */ 00314 if( dynamics.lgAdvection ) 00315 DynaStartZone(); 00316 00317 /* clear flag indicating the ionization convergence failures 00318 * have ever occurred in current zone 00319 conv.lgConvIonizThisZone = false; */ 00320 00321 /* this will say how many times the large H2 molecule has been called in this zone - 00322 * if not called (due to low H2 abundance) then not need to update its line arrays */ 00323 h2.nCallH2_this_zone = 0; 00324 00325 /* check whether zone thickness is too small relative to radius */ 00326 if( strcmp(dense.chDenseLaw,"GLOB") == 0 ) 00327 { 00328 ratio = radius.drad/(radius.glbdst + radius.drad); 00329 /* this only matters for globule command */ 00330 if( ratio < 1e-14 ) 00331 { 00332 radius.lgdR2Small = true; 00333 } 00334 else 00335 { 00336 radius.lgdR2Small = false; 00337 } 00338 } 00339 00340 00341 /* area factor, ratio of radius to out edge of this zone 00342 * relative to radius of illuminated face of cloud */ 00343 /*radius.r1r0sq = (float)POW2( 00344 (radius.Radius - radius.drad*radius.dRadSign/2.)/radius.rinner);*/ 00345 /*>>>chng 99 jan 25, definition of r1r0sq is now outer edge of current zone 00346 * relative to inner edge of cloud */ 00347 radius.r1r0sq = POW2( radius.Radius/radius.rinner ); 00348 00349 /* following only approximate, used for center of next zone */ 00350 radius.dRNeff = radius.drNext*geometry.FillFac; 00351 00352 /* this is the middle of the zone */ 00353 rad_middle_zone = radius.Radius - radius.drad/2.; 00354 00355 /* Radius is outer radius of this zone, so radius - drad is inner radius 00356 * rinner is inner radius of nebula; dVeff is the effective vol element 00357 * dVeff is vol rel to inner radius, so has units cm 00358 * for plane parallel dVeff is dReff */ 00359 if( radius.drad/radius.Radius > 0.01 ) 00360 { 00361 vin = ((radius.Radius - radius.drad)/radius.rinner)* 00362 (MIN2((radius.Radius-radius.drad),radius.CylindHigh)/radius.rinner)* 00363 (radius.Radius - radius.drad); 00364 00365 vout = (radius.Radius/radius.rinner)* 00366 (MIN2(radius.Radius,radius.CylindHigh)/radius.rinner)* 00367 radius.Radius; 00368 00369 /* default of iEmissPower is 2, observing the full sphere */ 00370 /* this is the usual case the full volume, just the difference in the two vols */ 00371 /* this will become the true vol when multiplied by the square of the inner radius */ 00372 radius.dVeff = (vout - vin)/3.*geometry.FillFac; 00373 } 00374 else 00375 { 00376 /* thin cell limit */ 00377 /* rad_middle_zone is the middle of the zone; 00378 * radius is always the OUTER edge of the current zone */ 00379 radius.dVeff = (rad_middle_zone/radius.rinner)*radius.drad*geometry.FillFac* 00380 (MIN2(rad_middle_zone,radius.CylindHigh)/radius.rinner); 00381 } 00382 00383 /* this is possible correction for slit projected onto resolved object - 00384 * this only happens when aperture command is entered. 00385 * NB not clear what happens when slit and cylinder are both used - 00386 * this would be very orientation specific */ 00387 00388 /* default of iEmissPower is 2, set to 0 is just aperture beam, 00389 * and to 1 if aperture slit set */ 00390 if( geometry.iEmissPower == 1 ) 00391 { 00392 /* this is long slit option, remove one power of radius */ 00393 radius.dVeff /= (rad_middle_zone/radius.rinner); 00394 } 00395 else if( geometry.iEmissPower == 0 ) 00396 { 00397 /* this is short slit option, remove two powers of radius */ 00398 radius.dVeff /= POW2(rad_middle_zone/radius.rinner); 00399 } 00400 00401 /* covering factor, default is unity */ 00402 radius.dVeff *= geometry.covgeo; 00403 00404 /* these are needed for line intensities in outward beam 00405 * lgSphere set, geometry.covrt usually 1, 0 when not lgSphere 00406 * so outward is 1 when lgSphere set 1/2 when not set, */ 00407 outwrd = (1. + geometry.covrt)/2.; 00408 00409 /*>>>chng 99 apr 23, from above to below */ 00410 /*radius.dVolOutwrd = outwrd*POW2( (radius.Radius-radius.drad_x_fillfac/2.)/radius.Radius) * 00411 radius.drad;*/ 00412 /* this includes covering fact, the effective vol,, and 1/r^2 dilution, 00413 * dReff includes filling factor. the covering factor part is 1 for sphere, 00414 * 1/2 for open */ 00415 /*radius.dVolOutwrd = outwrd*radius.Radius*radius.drad_x_fillfac/(radius.Radius + 00416 2.*radius.drad);*/ 00417 /* dReff from above, includes filling factor */ 00418 radius.dVolOutwrd = outwrd*radius.drad_x_fillfac; 00419 ASSERT( radius.dVolOutwrd > 0. ); 00420 00421 /* following will be used to "uncorrect" for this in lines when predicting continua 00422 radius.GeoDil = radius.Radius/(radius.Radius + 2.*radius.drad);*/ 00423 00424 /* this should multiply the line to become the net inward. geo part is 1/2 for 00425 * open geometry, 0 for closed. for case of isotropic emitter only this and dVolOutwrd 00426 * above are used */ 00427 radius.dVolReflec = (1. - outwrd)*radius.drad_x_fillfac*radius.r1r0sq; 00428 00429 if( geometry.lgSphere ) 00430 { 00431 /* inward beams do not go in when lgSphere set since geometry symmetric */ 00432 radius.BeamInIn = 0.; 00433 radius.BeamInOut = radius.Radius*radius.drad_x_fillfac/(radius.Radius + 00434 2.*radius.drad); 00435 } 00436 else 00437 { 00438 radius.BeamInIn = radius.drad_x_fillfac*radius.r1r0sq; 00439 00440 /* inward beams do not go out */ 00441 radius.BeamInOut = 0.; 00442 } 00443 /* factor for outwardly directed part of line */ 00444 radius.BeamOutOut = radius.Radius*radius.drad_x_fillfac/(radius.Radius + 00445 2.*radius.drad); 00446 00447 00448 DEBUG_EXIT( "ZoneStart()" ); 00449 return; 00450 } 00451 00452 void ZoneEnd(void) 00453 { 00454 long i; 00455 00456 DEBUG_ENTRY( "ZoneEnd()" ); 00457 00458 /*********************************************************************** 00459 * 00460 * correct rfield for attenuation from center of zone to inner edge 00461 * 00462 ***********************************************************************/ 00463 00464 /* radius is outer radius of this zone, this resets continuum to 00465 * flux at illuminated face of zone we have already computed. */ 00466 00467 /* opac.tmn defined in ZoneStart */ 00468 /* undilute radiation so that rfield is at face of zone */ 00469 /* NB, upper limit of sum includes [nflux] to test unit verification cell */ 00470 for( i=0; i <= rfield.nflux; i++ ) 00471 { 00472 rfield.flux[i] /= opac.tmn[i]; 00473 /* >>chng 03 nov 08, update SummedCon here since flux changed */ 00474 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i]; 00475 } 00476 00477 /* do advection if enabled */ 00478 if( dynamics.lgAdvection ) 00479 DynaEndZone(); 00480 00481 DEBUG_EXIT( "ZoneEnd()" ); 00482 return; 00483 }