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