cloudy trunk
|
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 /*tfidle update some temperature dependent variables */ 00004 /*tauff compute optical depth where cloud is thin to free-free and plasma freq */ 00005 /*velset set thermal velocities for all particles in gas */ 00006 #include "cddefines.h" 00007 #include "physconst.h" 00008 #include "opacity.h" 00009 #include "iso.h" 00010 #include "dense.h" 00011 #include "phycon.h" 00012 #include "stopcalc.h" 00013 #include "continuum.h" 00014 #include "trace.h" 00015 #include "rfield.h" 00016 #include "doppvel.h" 00017 #include "radius.h" 00018 #include "wind.h" 00019 #include "thermal.h" 00020 00021 /*tauff compute optical depth where cloud is thin to free-free and plasma freq */ 00022 STATIC void tauff(void); 00023 /*velset set thermal velocities for all particles in gas */ 00024 STATIC void velset(void); 00025 /* On first run, fill GauntFF with gaunt factors */ 00026 STATIC void FillGFF(void); 00027 /* Interpolate on GauntFF to calc gaunt at current temp, phycon.te */ 00028 STATIC realnum InterpolateGff( long charge , double ERyd ); 00029 STATIC int LinterpTable(realnum **t, realnum *v, long int lta, long int ltb, realnum x, realnum *a, long int *pipx); 00030 STATIC int LinterpVector(realnum **t, realnum *v, long lta , long ltb, realnum *yy , long ny, realnum **a); 00031 STATIC void fhunt(realnum *xx, long int n, realnum x, long int *j); 00032 00036 STATIC void tfidle( 00037 bool lgForceUpdate); 00038 00039 static long lgGffNotFilled = true; 00040 00041 #define N_TE_GFF 21L 00042 #define N_PHOTON_GFF 145L /* log(photon energy) = -8 to 10 in one-eighth dec steps */ 00043 static realnum ***GauntFF; 00044 static realnum **GauntFF_T; 00045 /* the array of logs of temperatures at which GauntFF is defined */ 00046 static realnum TeGFF[N_TE_GFF]; 00047 /* the array of logs of u at which GauntFF is defined */ 00048 static realnum PhoGFF[N_PHOTON_GFF]; 00049 00052 void TempChange( 00053 double TempNew , 00054 /* option to force update of all variables */ 00055 bool lgForceUpdate) 00056 { 00057 00058 DEBUG_ENTRY( "TempChange()" ); 00059 00060 /* set new temperature */ 00061 if( TempNew >= phycon.TEMP_LIMIT_HIGH ) 00062 { 00063 /* temp is too high */ 00064 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK," 00065 " is above the upper limit of the code, %.3eK.\n", 00066 TempNew , phycon.TEMP_LIMIT_HIGH ); 00067 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n"); 00068 00069 TempNew = phycon.TEMP_LIMIT_HIGH*0.999; 00070 lgAbort = true; 00071 } 00072 else if( TempNew <= phycon.TEMP_LIMIT_LOW ) 00073 { 00074 /* temp is too low */ 00075 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK," 00076 " is below the lower limit of the code, %.3eK.\n", 00077 TempNew , phycon.TEMP_LIMIT_LOW ); 00078 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n"); 00079 TempNew = phycon.TEMP_LIMIT_LOW*1.0001; 00080 lgAbort = true; 00081 } 00082 else if( TempNew < StopCalc.TeFloor ) 00083 { 00084 /* temperature floor option - 00085 * go to constant temperature calculation if temperature 00086 * falls below floor */ 00087 thermal.lgTemperatureConstant = true; 00088 thermal.ConstTemp = (realnum)StopCalc.TeFloor; 00089 phycon.te = thermal.ConstTemp; 00090 /*fprintf(ioQQQ,"DEBUG TempChange hit temp floor, setting const temp to %.3e\n", 00091 phycon.te );*/ 00092 } 00093 else 00094 { 00095 /* temp is within range */ 00096 phycon.te = TempNew; 00097 } 00098 00099 /* now update related variables */ 00100 tfidle( lgForceUpdate ); 00101 return; 00102 } 00106 void TempChange( 00107 double TempNew ) 00108 { 00109 00110 DEBUG_ENTRY( "TempChange()" ); 00111 00112 /* set new temperature */ 00113 if( TempNew >= phycon.TEMP_LIMIT_HIGH ) 00114 { 00115 /* temp is too high */ 00116 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK," 00117 " is above the upper limit of the code, %.3eK.\n", 00118 TempNew , phycon.TEMP_LIMIT_HIGH ); 00119 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n"); 00120 00121 TempNew = phycon.TEMP_LIMIT_HIGH*0.999; 00122 } 00123 else if( TempNew <= phycon.TEMP_LIMIT_LOW ) 00124 { 00125 /* temp is too low */ 00126 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK," 00127 " is below the lower limit of the code, %.3eK.\n", 00128 TempNew , phycon.TEMP_LIMIT_LOW ); 00129 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n"); 00130 TempNew = phycon.TEMP_LIMIT_LOW*1.0001; 00131 } 00132 else 00133 { 00134 /* temp is within range */ 00135 phycon.te = TempNew; 00136 } 00137 00138 /* now update related variables */ 00139 tfidle( false ); 00140 return; 00141 } 00142 00143 void tfidle( 00144 /* option to force update of all variables */ 00145 bool lgForceUpdate) 00146 { 00147 static double tgffused=-1., 00148 tgffsued2=-1.; 00149 static long int nff_defined=-1; 00150 static long maxion = 0, oldmaxion = 0; 00151 static double ttused = 0.; 00152 static double edused = 0.; 00153 static bool lgZLogSet = false; 00154 bool lgGauntF; 00155 long int ion; 00156 long int i, 00157 nelem, 00158 if1, 00159 ipTe, 00160 ret; 00161 double fac, 00162 fanu; 00163 00164 DEBUG_ENTRY( "tfidle()" ); 00165 00166 /* called with lgForceUpdate true in zero.c, when we must update everything */ 00167 if( lgForceUpdate ) 00168 { 00169 ttused = -1.; 00170 edused = 0.; 00171 } 00172 00173 /* check that eden not negative */ 00174 if( dense.eden <= 0. ) 00175 { 00176 fprintf( ioQQQ, "tfidle called with a zero or negative electron density,%10.2e\n", 00177 dense.eden ); 00178 TotalInsanity(); 00179 } 00180 00181 /* check that temperature not negative */ 00182 if( phycon.te <= 0. ) 00183 { 00184 fprintf( ioQQQ, "tfidle called with a negative electron temperature,%10.2e\n", 00185 phycon.te ); 00186 TotalInsanity(); 00187 } 00188 00189 /* only time only, set up array of logs of charge squared */ 00190 if( !lgZLogSet ) 00191 { 00192 for( nelem=0; nelem<LIMELM; ++nelem ) 00193 { 00194 /* this array is used to modify the log temperature array 00195 * defined below, for hydrogenic species of charge nelem+1 */ 00196 phycon.sqlogz[nelem] = log10( POW2(nelem+1.) ); 00197 } 00198 lgZLogSet = true; 00199 } 00200 00201 if( ! fp_equal( phycon.te, ttused ) ) 00202 { 00203 ttused = phycon.te; 00204 thermal.te_update = phycon.te; 00205 /* current temperature in various units */ 00206 phycon.te_eV = phycon.te/EVDEGK; 00207 phycon.te_ryd = phycon.te/TE1RYD; 00208 phycon.te_wn = phycon.te / T1CM; 00209 00210 phycon.tesqrd = POW2(phycon.te); 00211 phycon.sqrte = sqrt(phycon.te); 00212 thermal.halfte = 0.5/phycon.te; 00213 thermal.tsq1 = 1./phycon.tesqrd; 00214 phycon.te32 = phycon.te*phycon.sqrte; 00215 phycon.teinv = 1./phycon.te; 00216 00217 phycon.alogte = log10(phycon.te); 00218 phycon.alnte = log(phycon.te); 00219 00220 phycon.telogn[0] = phycon.alogte; 00221 for( i=1; i < 7; i++ ) 00222 { 00223 phycon.telogn[i] = phycon.telogn[i-1]*phycon.telogn[0]; 00224 } 00225 00226 phycon.te10 = pow(phycon.te,0.10); 00227 phycon.te20 = phycon.te10 * phycon.te10; 00228 phycon.te30 = phycon.te20 * phycon.te10; 00229 phycon.te40 = phycon.te30 * phycon.te10; 00230 phycon.te70 = phycon.sqrte * phycon.te20; 00231 phycon.te90 = phycon.te70 * phycon.te20; 00232 00233 phycon.te01 = pow(phycon.te,0.01); 00234 phycon.te02 = phycon.te01 * phycon.te01; 00235 phycon.te03 = phycon.te02 * phycon.te01; 00236 phycon.te04 = phycon.te02 * phycon.te02; 00237 phycon.te05 = phycon.te03 * phycon.te02; 00238 phycon.te07 = phycon.te05 * phycon.te02; 00239 00240 phycon.te001 = pow(phycon.te,0.001); 00241 phycon.te002 = phycon.te001 * phycon.te001; 00242 phycon.te003 = phycon.te002 * phycon.te001; 00243 phycon.te004 = phycon.te002 * phycon.te002; 00244 phycon.te005 = phycon.te003 * phycon.te002; 00245 phycon.te007 = phycon.te005 * phycon.te002; 00246 /*>>>chng 06 June 30 -Humeshkar Nemala*/ 00247 phycon.te0001 = pow(phycon.te ,0.0001); 00248 phycon.te0002 = phycon.te0001 * phycon.te0001; 00249 phycon.te0003 = phycon.te0002 * phycon.te0001; 00250 phycon.te0004 = phycon.te0002 * phycon.te0002; 00251 phycon.te0005 = phycon.te0003 * phycon.te0002; 00252 phycon.te0007 = phycon.te0005 * phycon.te0002; 00253 00254 } 00255 00256 /* >>>chng 99 nov 23, removed this line, so back to old method of h coll */ 00257 /* used for hydrogenic collisions */ 00258 /* 00259 * following electron density has approximate correction for neutrals 00260 * corr of hi*1.7e-4 accounts for col ion by HI; 00261 * >>refer H0 correction for collisional contribution Drawin, H.W. 1969, Zs Phys 225, 483. 00262 * also quoted in Dalgarno & McCray 1972 00263 * extensive discussion of this in 00264 *>>refer H0 collisions Lambert, D.L. 00265 * used EdenHCorr instead 00266 * edhi = eden + hi * 1.7e-4 00267 */ 00268 dense.EdenHCorr = dense.eden + 00269 /* dense.HCorrFac is unity by default and changed with the set HCOR command */ 00270 dense.xIonDense[ipHYDROGEN][0]*1.7e-4 * dense.HCorrFac; 00271 00272 /* this is parallel version for H0 onto H0 collisions - for now the same, but 00273 * this may not be correct */ 00274 dense.EdenHontoHCorr = dense.EdenHCorr; 00275 00276 /*>>chng 93 jun 04, 00277 * term with hi added June 4, 93, to account for warm pdr */ 00278 /* >>chng 05 jan 05, Will Henney noticed that 1.e-4 used here is not same as 00279 * 1.7e-4 used for EdenHCorr, which had rewritten the expression. 00280 * change so that edensqte uses EdenHCorr rather than reevaluating */ 00281 /*dense.edensqte = ((dense.eden + dense.xIonDense[ipHYDROGEN][0]/1e4)/phycon.sqrte);*/ 00282 dense.edensqte = dense.EdenHCorr/phycon.sqrte; 00283 dense.cdsqte = dense.edensqte*COLL_CONST; 00284 00285 if( fabs(1.-edused/phycon.te)>=0.00001 ) 00286 { 00287 edused = dense.eden; 00288 dense.SqrtEden = sqrt(dense.eden); 00289 } 00290 00291 /* finally reset velocities */ 00292 /* find line widths for thermal and turbulent motion 00293 * CLOUDY uses line center optical depths, =0.015 F / DELTA NU */ 00294 velset(); 00295 00296 /* rest have to do with radiation field which may not be defined yet */ 00297 if( !lgRfieldMalloced ) 00298 { 00299 return; 00300 } 00301 00302 /* correction factors for induced recombination, 00303 * also used as Boltzmann factors 00304 * check for zero is because ContBoltz is zeroed out in initialization 00305 * of code, its possible this is a constant density grid of models 00306 * in which the code is called as a subroutine */ 00307 /* >>chng 01 aug 21, must also test on size of continuum nflux because 00308 * conintitemp can increase nflux then call this routine, although 00309 * temp may not have changed */ 00310 if( fabs(1.-tgffused/phycon.te)>=0.00001 /* tgffused != phycon.te */ 00311 || rfield.ContBoltz[0] <= 0. || nff_defined<rfield.nflux ) 00312 { 00313 tgffused = phycon.te; 00314 fac = TE1RYD/phycon.te; 00315 i = 0; 00316 fanu = fac*rfield.anu[i]; 00317 /* SEXP_LIMIT is 84 in cddefines.h - it is the -ln of smallest number 00318 * that sexp can handle, and is used elsewhere in the code. 00319 * atom_level2 uses ContBoltz to see whether the rates will be significant. 00320 * If the numbers did not agree then this test would be flawed, resulting in 00321 * div by zero */ 00322 while( i < rfield.nupper && fanu < SEXP_LIMIT ) 00323 { 00324 rfield.ContBoltz[i] = exp(-fanu); 00325 ++i; 00326 /* this is Boltzmann factor for NEXT cell */ 00327 fanu = fac*rfield.anu[i]; 00328 } 00329 /* ipMaxBolt is number of cells defined, so defined up through ipMaxBolt-1 */ 00330 rfield.ipMaxBolt = i; 00331 00332 /* zero out remainder */ 00333 /* >>chng 01 apr 14, upper limit has been ipMaxBolt+1, off by one */ 00334 for( i=rfield.ipMaxBolt; i < rfield.nupper; i++ ) 00335 { 00336 rfield.ContBoltz[i] = 0.; 00337 } 00338 } 00339 00340 /* find frequency where thin to bremsstrahlung or plasma frequency */ 00341 tauff(); 00342 00343 oldmaxion = maxion; 00344 maxion = 0; 00345 00346 /* Find highest maximum stage of ionization */ 00347 for( nelem = 0; nelem < LIMELM; nelem++ ) 00348 { 00349 maxion = MAX2( maxion, dense.IonHigh[nelem] ); 00350 } 00351 00352 /* reevaluate if temperature or number of cells has changed */ 00353 /* >>chng 03 sep 26, following had been 0.1, causing jumps when evaluated, 00354 * changed to 0.02 */ 00355 #define LIM 0.02 00356 if( fabs(1.-tgffsued2/phycon.te) > LIM || 00357 /* this test - reevaluate if upper bound of defined values is 00358 * above nflux, the highest point. This only triggers in 00359 * large grids when continuum source gets harder */ 00360 nff_defined<rfield.nflux 00361 /* this occurs when now have more stages of ionization than in previous defn */ 00362 || maxion>oldmaxion ) 00363 { 00364 static bool lgFirstRunDone = false; 00365 long lowion; 00366 /* >>chng 02 jan 10, only reevaluate needed states of ionization */ 00367 if( fabs(1.-tgffsued2/phycon.te) <= LIM && nff_defined>=rfield.nflux && 00368 maxion>oldmaxion ) 00369 { 00370 /* this case temperature did not change by much, but upper 00371 * stage of ionization increased. only need to evaluate 00372 * stages that have not been already done */ 00373 lowion = oldmaxion; 00374 } 00375 else 00376 { 00377 /* temperature changed - do them all */ 00378 lowion = 1; 00379 } 00380 00381 /* if1 will certainly be set to some positive value in gffsub */ 00382 if1 = 1; 00383 00384 /* chng 02 may 16, by Ryan...one gaunt factor array for all charges */ 00385 /* First index is EFFECTIVE CHARGE! */ 00386 /* highest stage of ionization is LIMELM, 00387 * index goes from 1 to LIMELM */ 00388 00389 nff_defined = rfield.nflux; 00390 00391 ASSERT( if1 >= 0 ); 00392 00393 tgffsued2 = phycon.te; 00394 lgGauntF = true; 00395 00396 /* new gaunt factors */ 00397 if( lgGffNotFilled ) 00398 { 00399 FillGFF(); 00400 } 00401 00402 if( lgFirstRunDone == false ) 00403 { 00404 maxion = LIMELM; 00405 lgFirstRunDone = true; 00406 } 00407 00408 /* >> chng 03 jan 23, rjrw -- move a couple of loops down into 00409 * subroutines, and make those subroutines generic 00410 * (i.e. dependences only on arguments, may be useful elsewhere...) */ 00411 00412 /* Make Gaunt table for new temperature */ 00413 ipTe = -1; 00414 for( ion=1; ion<=LIMELM; ion++ ) 00415 { 00416 /* Interpolate for all tabulated photon energies at this temperature */ 00417 ret = LinterpTable(GauntFF[ion], TeGFF, N_PHOTON_GFF, N_TE_GFF, (realnum)phycon.alogte, GauntFF_T[ion], &ipTe); 00418 if(ret == -1) 00419 { 00420 fprintf(ioQQQ," LinterpTable for GffTable called with te out of bounds \n"); 00421 cdEXIT(EXIT_FAILURE); 00422 } 00423 } 00424 00425 /* Interpolate for all ions at required photon energies -- similar 00426 * to LinterpTable, but not quite similar enough... */ 00427 /* >>chng 04 jun 30, change nflux to nupper */ 00428 ret = LinterpVector(GauntFF_T+lowion, PhoGFF, maxion-lowion+1, N_PHOTON_GFF, 00429 rfield.anulog, rfield.nupper,/*rfield.nflux,*/ rfield.gff + lowion); 00430 if(ret == -1) 00431 { 00432 fprintf(ioQQQ," LinterpVector for GffTable called with photon energy out of bounds \n"); 00433 cdEXIT(EXIT_FAILURE); 00434 } 00435 } 00436 else 00437 { 00438 /* this is flag that would have been set in gffsub, and 00439 * printed in debug statement below. We are not evaluating 00440 * so set to -1 */ 00441 if1 = -1; 00442 lgGauntF = false; 00443 } 00444 00445 if( trace.lgTrace && trace.lgTrGant ) 00446 { 00447 fprintf( ioQQQ, " tfidle; gaunt factors?" ); 00448 fprintf( ioQQQ, "%2c", TorF(lgGauntF) ); 00449 00450 fprintf( ioQQQ, "%2f g 1 2=%10.2e%9.2ld flag%2f guv(1,n)%10.2e\n", 00451 rfield.gff[1][0], rfield.gff[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1], 00452 if1, rfield.gff[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]], 00453 rfield.gff[1][rfield.nflux-1] ); 00454 } 00455 return; 00456 } 00457 00458 /*tauff compute optical depth where cloud is thin to free-free and plasma freq */ 00459 STATIC void tauff(void) 00460 { 00461 realnum fac; 00462 00463 /* simply return if space not yet allocated */ 00464 if( !lgOpacMalloced ) 00465 return; 00466 00467 DEBUG_ENTRY( "tauff()" ); 00468 00469 /* routine sets variable ipEnergyBremsThin, index for lowest cont cell that is optically thin */ 00470 /* find frequency where continuum thin to free-free */ 00471 while( rfield.ipEnergyBremsThin < rfield.nflux && 00472 opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] >= 1. ) 00473 { 00474 ++rfield.ipEnergyBremsThin; 00475 } 00476 00477 /* TFF will be frequency where cloud becomes optically thin to bremsstrahlung 00478 * >>chng 96 may 7, had been 2, change as per Kevin Volk bug report */ 00479 if( rfield.ipEnergyBremsThin > 1 && opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] > 0.001 ) 00480 { 00481 /* tau can be zero when plasma frequency is within energy grid, */ 00482 fac = (1.f - opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1])/(opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] - 00483 opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1]); 00484 fac = MAX2(fac,0.f); 00485 rfield.EnergyBremsThin = rfield.anu[rfield.ipEnergyBremsThin-1] + rfield.widflx[rfield.ipEnergyBremsThin-1]*fac; 00486 } 00487 else 00488 { 00489 rfield.EnergyBremsThin = 0.f; 00490 } 00491 00492 /* now evaluate the plasma frequency */ 00493 rfield.plsfrq = (realnum)(2.729e-12*sqrt(dense.eden*1.2)); 00494 00495 if( rfield.ipPlasma > 0 ) 00496 { 00497 /* >>chng 02 jul 25, increase index for plasma frequency until within proper cell */ 00498 while( rfield.plsfrq > rfield.anu[rfield.ipPlasma]+rfield.widflx[rfield.ipPlasma]/2. ) 00499 ++rfield.ipPlasma; 00500 00501 /* >>chng 02 jul 25, decrease index for plasma frequency until within proper cell */ 00502 while( rfield.ipPlasma>2 && rfield.plsfrq < rfield.anu[rfield.ipPlasma]-rfield.widflx[rfield.ipPlasma]/2. ) 00503 --rfield.ipPlasma; 00504 } 00505 00506 /* also remember the largest plasma frequency we encounter */ 00507 rfield.plsfrqmax = MAX2(rfield.plsfrqmax, rfield.plsfrq); 00508 00509 /* is plasma frequency within energy grid? */ 00510 if( rfield.plsfrq > rfield.anu[0] ) 00511 { 00512 rfield.lgPlasNu = true; 00513 } 00514 00515 /* >>chng 96 jul 15, did not include plasma freq before 00516 * function returns larger of these two frequencies */ 00517 rfield.EnergyBremsThin = MAX2(rfield.plsfrq,rfield.EnergyBremsThin); 00518 00519 /* now increment ipEnergyBremsThin still further, until above plasma frequency */ 00520 while( rfield.ipEnergyBremsThin < rfield.nflux && 00521 rfield.anu[rfield.ipEnergyBremsThin] <= rfield.EnergyBremsThin ) 00522 { 00523 ++rfield.ipEnergyBremsThin; 00524 } 00525 return; 00526 } 00527 00528 /*velset set thermal velocities for all particles in gas */ 00529 STATIC void velset(void) 00530 { 00531 long int nelem; 00532 double turb2; 00533 00534 DEBUG_ENTRY( "velset()" ); 00535 00536 /* usually TurbVel =0, reset with turbulence command 00537 * cm/s here, but was entered in km/s with command */ 00538 turb2 = POW2(DoppVel.TurbVel); 00539 00540 /* this is option to dissipate the turbulence. DispScale is entered with 00541 * dissipate keyword on turbulence command. The velocity is reduced here, 00542 * by an assumed exponential scale, and also adds heat */ 00543 if( DoppVel.DispScale > 0. ) 00544 { 00545 /* square of exp depth dependence */ 00546 turb2 *= sexp( 2.*radius.depth / DoppVel.DispScale ); 00547 } 00548 00549 /* in case of D-Critical flow include initial velocity as 00550 * a component of turbulence */ 00551 if( wind.windv0 < 0. ) 00552 { 00553 turb2 += POW2(wind.windv0); 00554 } 00555 00556 /* computes one doppler width, in cm/sec, 00557 * for each element with atomic number the array index*/ 00558 for( nelem=0; nelem < LIMELM; nelem++ ) 00559 { 00560 DoppVel.doppler[nelem] = 00561 /*(realnum)sqrt(1.651e8*phycon.te/dense.AtomicWeight[nelem]+*/ 00562 /* >>chng 00 may 02 to physical constants */ 00563 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/dense.AtomicWeight[nelem]+ 00564 turb2); 00565 /* this is average (NOT rms) particle speed for Maxwell distribution, Mihalas 70, 9-70 */ 00566 DoppVel.AveVel[nelem] = sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT*phycon.te/dense.AtomicWeight[nelem]); 00567 } 00568 00569 /* DoppVel.doppler[LIMELM] is CO, vector is dim LIMELM+1 */ 00570 /* C12O16 */ 00571 DoppVel.doppler[LIMELM] = 00572 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/ 00573 (dense.AtomicWeight[5]+dense.AtomicWeight[7]) + turb2); 00574 DoppVel.AveVel[LIMELM] = 00575 sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/ 00576 (dense.AtomicWeight[5]+dense.AtomicWeight[7]) + turb2); 00577 00578 /* C13O16 */ 00579 DoppVel.doppler[LIMELM+1] = 00580 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/ 00581 (dense.AtomicWeight[5]*13./12.+dense.AtomicWeight[7]) + turb2); 00582 DoppVel.AveVel[LIMELM+1] = 00583 sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/ 00584 (dense.AtomicWeight[5]*13./12.+dense.AtomicWeight[7]) + turb2); 00585 00586 /* H2 */ 00587 DoppVel.doppler[LIMELM+2] = 00588 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/ 00589 (2.*dense.AtomicWeight[0]) + turb2); 00590 DoppVel.AveVel[LIMELM+2] = 00591 sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/ 00592 (2.*dense.AtomicWeight[0]) + turb2); 00593 return; 00594 } 00595 00596 STATIC void FillGFF( void ) 00597 { 00598 00599 long i,i1,i2,i3,j,charge,GffMAGIC = 80321; 00600 double Temp, photon; 00601 bool lgEOL; 00602 00603 DEBUG_ENTRY( "FillGFF()" ); 00604 00605 # define chLine_LENGTH 1000 00606 char chLine[chLine_LENGTH]; 00607 00608 FILE *ioDATA; 00609 00610 for( i = 0; i < N_TE_GFF; i++ ) 00611 { 00612 TeGFF[i] = 0.5f*i; 00613 } 00614 /* >>chng 06 feb 14, assert thrown at T == 1e10 K, Ryan Porter proposes this fix */ 00615 TeGFF[N_TE_GFF-1] += 0.01f; 00616 00617 for( i = 0; i< N_PHOTON_GFF; i++ ) 00618 { 00619 PhoGFF[i] = 0.125f*i - 8.f; 00620 } 00621 00622 GauntFF = (realnum***)MALLOC( sizeof(realnum**)*(unsigned)(LIMELM+2) ); 00623 for( i = 1; i <= LIMELM; i++ ) 00624 { 00625 GauntFF[i] = (realnum**)MALLOC( sizeof(realnum*)*(unsigned)N_PHOTON_GFF ); 00626 for( j = 0; j < N_PHOTON_GFF; j++ ) 00627 { 00628 GauntFF[i][j] = (realnum*)MALLOC( sizeof(realnum)*(unsigned)N_TE_GFF ); 00629 } 00630 } 00631 00632 GauntFF_T = (realnum**)MALLOC( sizeof(realnum*)*(unsigned)(LIMELM+2) ); 00633 for( i = 1; i <= LIMELM; i++ ) 00634 { 00635 GauntFF_T[i] = (realnum*)MALLOC( sizeof(realnum)*(unsigned)N_PHOTON_GFF ); 00636 } 00637 00638 if( !rfield.lgCompileGauntFF ) 00639 { 00640 if( trace.lgTrace ) 00641 fprintf( ioQQQ," FillGFF opening gauntff.dat:"); 00642 00643 try 00644 { 00645 ioDATA = open_data( "gauntff.dat", "r" ); 00646 } 00647 catch( cloudy_exit ) 00648 { 00649 fprintf( ioQQQ, " Defaulting to on-the-fly computation, "); 00650 fprintf( ioQQQ, "but the code runs much faster if you compile gauntff.dat!\n"); 00651 ioDATA = NULL; 00652 } 00653 00654 if( ioDATA == NULL ) 00655 { 00656 /* Do on the fly computation of Gff's instead. */ 00657 for( charge=1; charge<=LIMELM; charge++ ) 00658 { 00659 for( i=0; i<N_PHOTON_GFF; i++ ) 00660 { 00661 photon = pow((realnum)10.f,PhoGFF[i]); 00662 for(j=0; j<N_TE_GFF; j++) 00663 { 00664 Temp = pow((realnum)10.f,TeGFF[j]); 00665 GauntFF[charge][i][j] = (realnum)cont_gaunt_calc( Temp, (double)charge, photon ); 00666 } 00667 } 00668 } 00669 } 00670 else 00671 { 00672 /* check that magic number is ok */ 00673 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00674 { 00675 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n"); 00676 cdEXIT(EXIT_FAILURE); 00677 } 00678 i = 1; 00679 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00680 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00681 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00682 00683 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF ) 00684 { 00685 fprintf( ioQQQ, 00686 " FillGFF: the version of gauntff.dat is not the current version.\n" ); 00687 fprintf( ioQQQ, 00688 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" , 00689 GffMAGIC , 00690 N_PHOTON_GFF, 00691 N_TE_GFF, 00692 i1 , i2 , i3 ); 00693 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00694 fprintf( ioQQQ, 00695 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" ); 00696 cdEXIT(EXIT_FAILURE); 00697 } 00698 00699 /* now read in the data */ 00700 for( charge = 1; charge <= LIMELM; charge++ ) 00701 { 00702 for( i = 0; i<N_PHOTON_GFF; i++ ) 00703 { 00704 /* get next line image */ 00705 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00706 { 00707 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n"); 00708 cdEXIT(EXIT_FAILURE); 00709 } 00710 /* each line starts with charge and energy level ( index in rfield ) */ 00711 i3 = 1; 00712 i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL); 00713 i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL); 00714 /* check that these numbers are correct */ 00715 if( i1!=charge || i2!=i ) 00716 { 00717 fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n"); 00718 fprintf( ioQQQ, 00719 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" ); 00720 cdEXIT(EXIT_FAILURE); 00721 } 00722 00723 /* loop over temperatures to produce array of free free gaunt factors */ 00724 for(j = 0; j < N_TE_GFF; j++) 00725 { 00726 GauntFF[charge][i][j] = (realnum)FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL); 00727 00728 if( lgEOL ) 00729 { 00730 fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n"); 00731 fprintf( ioQQQ, 00732 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" ); 00733 cdEXIT(EXIT_FAILURE); 00734 } 00735 } 00736 } 00737 00738 } 00739 00740 /* check that magic number is ok */ 00741 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00742 { 00743 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n"); 00744 cdEXIT(EXIT_FAILURE); 00745 } 00746 i = 1; 00747 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00748 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00749 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00750 00751 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF ) 00752 { 00753 fprintf( ioQQQ, 00754 " FillGFF: the version of gauntff.dat is not the current version.\n" ); 00755 fprintf( ioQQQ, 00756 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" , 00757 GffMAGIC , 00758 N_PHOTON_GFF, 00759 N_TE_GFF, 00760 i1 , i2 , i3 ); 00761 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00762 fprintf( ioQQQ, 00763 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" ); 00764 cdEXIT(EXIT_FAILURE); 00765 } 00766 00767 /* close the data file */ 00768 fclose( ioDATA ); 00769 } 00770 if( trace.lgTrace ) 00771 fprintf( ioQQQ," - opened and read ok.\n"); 00772 } 00773 else 00774 { 00775 /* option to create table of gaunt factors, 00776 * executed with the compile gaunt command */ 00777 FILE *ioGFF; 00778 00779 /*GffMAGIC the magic number for the table of recombination coefficients */ 00780 ioGFF = open_data( "gauntff.dat" , "w", AS_LOCAL_ONLY ); 00781 fprintf(ioGFF,"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n", 00782 GffMAGIC , 00783 N_PHOTON_GFF, 00784 N_TE_GFF, 00785 N_PHOTON_GFF, 00786 N_TE_GFF ); 00787 00788 for( charge = 1; charge <= LIMELM; charge++ ) 00789 { 00790 for( i=0; i < N_PHOTON_GFF; i++ ) 00791 { 00792 fprintf(ioGFF, "%li\t%li", charge, i ); 00793 /* loop over temperatures to produce array of gaunt factors */ 00794 for(j = 0; j < N_TE_GFF; j++) 00795 { 00796 /* Store gaunt factor in N_TE_GFF half dec steps */ 00797 Temp = pow((realnum)10.f,TeGFF[j]); 00798 photon = pow((realnum)10.f,PhoGFF[i]); 00799 GauntFF[charge][i][j] = (realnum)cont_gaunt_calc( Temp, (double)charge, photon ); 00800 fprintf(ioGFF, "\t%f", GauntFF[charge][i][j] ); 00801 } 00802 fprintf(ioGFF, "\n" ); 00803 } 00804 } 00805 00806 /* end the file with the same information */ 00807 fprintf(ioGFF,"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n", 00808 GffMAGIC , 00809 N_PHOTON_GFF, 00810 N_TE_GFF, 00811 N_PHOTON_GFF, 00812 N_TE_GFF ); 00813 00814 fclose( ioGFF ); 00815 00816 fprintf( ioQQQ, "FillGFF: compilation complete, gauntff.dat created.\n" ); 00817 } 00818 00819 lgGffNotFilled = false; 00820 00821 { 00822 /* We have already checked the validity of cont_gaunt_calc in sanitycheck.c. 00823 * Now we check to see if the InterpolateGff routine also works correctly. */ 00824 /*@-redef@*/ 00825 enum {DEBUG_LOC=false}; 00826 /* if set true there will be two problems at very low temps */ 00827 /*@+redef@*/ 00828 if( DEBUG_LOC ) 00829 { 00830 double gaunt, error; 00831 double tempsave = phycon.te; 00832 long logu, loggamma2; 00833 00834 for( logu=-4; logu<=4; logu++) 00835 { 00836 /* Uncommenting each of the three print statements in this bit 00837 * will produce a nice table comparable to Table 2 of 00838 * >>refer free-free gaunts Sutherland, R.S., 1998, MNRAS, 300, 321-330 */ 00839 /* fprintf(ioQQQ,"%li\t", logu);*/ 00840 for(loggamma2=-4; loggamma2<=4; loggamma2++) 00841 { 00842 double SutherlandGff[9][9]= 00843 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008}, 00844 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041}, 00845 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771}, 00846 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747}, 00847 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237}, 00848 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202}, 00849 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355}, 00850 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680}, 00851 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}}; 00852 00853 phycon.te = (TE1RYD/pow(10.,(double)loggamma2)); 00854 phycon.alogte = log10(phycon.te); 00855 gaunt = InterpolateGff( 1, pow(10.,(double)(logu-loggamma2)) ); 00856 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt; 00857 /*fprintf(ioQQQ,"%1.3f\t", gaunt);*/ 00858 if( error>0.05 ) 00859 { 00860 fprintf(ioQQQ," tfidle found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n", 00861 logu, loggamma2, error ); 00862 } 00863 } 00864 /*fprintf(ioQQQ,"\n");*/ 00865 } 00866 phycon.te = tempsave; 00867 phycon.alogte = log10(phycon.te); 00868 } 00869 } 00870 00871 return; 00872 } 00873 00874 /* Interpolate Gff at some temperature */ 00875 STATIC realnum InterpolateGff( long charge , double ERyd ) 00876 { 00877 double GauntAtLowerPho, GauntAtUpperPho; 00878 static long int ipTe=-1, ipPho=-1; 00879 double gaunt = 0.; 00880 double xmin , xmax; 00881 long i; 00882 00883 DEBUG_ENTRY( "InterpolateGff()" ); 00884 00885 if( ipTe < 0 ) 00886 { 00887 /* te totally unknown */ 00888 if( phycon.alogte < TeGFF[0] || phycon.alogte > TeGFF[N_TE_GFF-1] ) 00889 { 00890 fprintf(ioQQQ," InterpolateGff called with te out of bounds \n"); 00891 cdEXIT(EXIT_FAILURE); 00892 } 00893 /* now search for temperature */ 00894 for( i=0; i<N_TE_GFF-1; ++i ) 00895 { 00896 if( phycon.alogte > TeGFF[i] && phycon.alogte <= TeGFF[i+1] ) 00897 { 00898 /* found the temperature - use it */ 00899 ipTe = i; 00900 break; 00901 } 00902 } 00903 ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1) ); 00904 00905 } 00906 else if( phycon.alogte < TeGFF[ipTe] ) 00907 { 00908 /* temp is too low, must also lower ipTe */ 00909 ASSERT( phycon.alogte > TeGFF[0] ); 00910 /* decrement ipTe until it is correct */ 00911 while( phycon.alogte < TeGFF[ipTe] && ipTe > 0) 00912 --ipTe; 00913 } 00914 else if( phycon.alogte > TeGFF[ipTe+1] ) 00915 { 00916 /* temp is too high */ 00917 ASSERT( phycon.alogte <= TeGFF[N_TE_GFF-1] ); 00918 /* increment ipTe until it is correct */ 00919 while( phycon.alogte > TeGFF[ipTe+1] && ipTe < N_TE_GFF-1) 00920 ++ipTe; 00921 } 00922 00923 ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1) ); 00924 00925 /* ipTe should now be valid */ 00926 ASSERT( phycon.alogte >= TeGFF[ipTe] && phycon.alogte <= TeGFF[ipTe+1] && ipTe < N_TE_GFF-1 ); 00927 00928 /***************/ 00929 /* This bit is completely analogous to the above, but for the photon vector instead of temp. */ 00930 if( ipPho < 0 ) 00931 { 00932 if( log10(ERyd) < PhoGFF[0] || log10(ERyd) > PhoGFF[N_PHOTON_GFF-1] ) 00933 { 00934 fprintf(ioQQQ," InterpolateGff called with photon energy out of bounds \n"); 00935 cdEXIT(EXIT_FAILURE); 00936 } 00937 for( i=0; i<N_PHOTON_GFF-1; ++i ) 00938 { 00939 if( log10(ERyd) > PhoGFF[i] && log10(ERyd) <= PhoGFF[i+1] ) 00940 { 00941 ipPho = i; 00942 break; 00943 } 00944 } 00945 ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1) ); 00946 00947 } 00948 else if( log10(ERyd) < PhoGFF[ipPho] ) 00949 { 00950 ASSERT( log10(ERyd) >= PhoGFF[0] ); 00951 while( log10(ERyd) < PhoGFF[ipPho] && ipPho > 0) 00952 --ipPho; 00953 } 00954 else if( log10(ERyd) > PhoGFF[ipPho+1] ) 00955 { 00956 ASSERT( log10(ERyd) <= PhoGFF[N_PHOTON_GFF-1] ); 00957 while( log10(ERyd) > PhoGFF[ipPho+1] && ipPho < N_PHOTON_GFF-1) 00958 ++ipPho; 00959 } 00960 ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1) ); 00961 ASSERT( log10(ERyd)>=PhoGFF[ipPho] 00962 && log10(ERyd)<=PhoGFF[ipPho+1] && ipPho<N_PHOTON_GFF-1 ); 00963 00964 /* Calculate the answer...must interpolate on two variables. 00965 * First interpolate on T, at both the lower and upper photon energies. 00966 * Then interpolate between these results for the right photon energy. */ 00967 00968 GauntAtLowerPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) * 00969 (GauntFF[charge][ipPho][ipTe+1] - GauntFF[charge][ipPho][ipTe]) + GauntFF[charge][ipPho][ipTe]; 00970 00971 GauntAtUpperPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) * 00972 (GauntFF[charge][ipPho+1][ipTe+1] - GauntFF[charge][ipPho+1][ipTe]) + GauntFF[charge][ipPho+1][ipTe]; 00973 00974 gaunt = (log10(ERyd) - PhoGFF[ipPho]) / (PhoGFF[ipPho+1] - PhoGFF[ipPho]) * 00975 (GauntAtUpperPho - GauntAtLowerPho) + GauntAtLowerPho; 00976 00977 xmax = MAX4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1], 00978 GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] ); 00979 ASSERT( gaunt <= xmax ); 00980 00981 xmin = MIN4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1], 00982 GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] ); 00983 ASSERT( gaunt >= xmin ); 00984 00985 ASSERT( gaunt > 0. ); 00986 00987 return (realnum)gaunt; 00988 } 00989 00990 /* Interpolate in table t[lta][ltb], with physical values for the 00991 second index given by v[ltb], for values x, and put results in 00992 a[lta]; store the index found if that's useful; assumes v[] is 00993 sorted */ 00994 STATIC int LinterpTable(realnum **t, realnum *v, long int lta, long int ltb, realnum x, realnum *a, long int *pipx) 00995 { 00996 long int ipx=-1; 00997 realnum frac; 00998 long i; 00999 01000 DEBUG_ENTRY( "LinterpTable()" ); 01001 01002 if(pipx != NULL) 01003 ipx = *pipx; 01004 01005 fhunt (v,ltb,x,&ipx); /* search for index */ 01006 if(pipx != NULL) 01007 *pipx = ipx; 01008 01009 if( ipx == -1 || ipx == ltb ) 01010 { 01011 return -1; 01012 } 01013 01014 ASSERT( (ipx >=0) && (ipx < ltb-1) ); 01015 ASSERT( x >= v[ipx] && x <= v[ipx+1]); 01016 01017 frac = (x - v[ipx]) / (v[ipx+1] - v[ipx]); 01018 for( i=0; i<lta; i++ ) 01019 { 01020 a[i] = frac*t[i][ipx+1]+(1.f-frac)*t[i][ipx]; 01021 } 01022 01023 return 0; 01024 } 01025 01026 /* Interpolate in table t[lta][ltb], with physical values for the second index given by v[ltb], 01027 for values yy[ny], and put results in a[lta][ly] */ 01028 STATIC int LinterpVector(realnum **t, realnum *v, long lta , long ltb, realnum *yy, long ly, realnum **a) 01029 { 01030 realnum yl, frac; 01031 long i, j, n; 01032 01033 DEBUG_ENTRY( "LinterpVector()" ); 01034 01035 if( yy[0] < v[0] || yy[ly-1] > v[ltb-1] ) 01036 { 01037 return -1; 01038 } 01039 01040 n = 0; 01041 yl = yy[n]; 01042 for(j = 1; j < ltb && n < ly; j++) { 01043 while (yl < v[j] && n < ly) { 01044 frac = (yl-v[j-1])/(v[j]-v[j-1]); 01045 for(i = 0; i < lta; i++) 01046 a[i][n] = frac*t[i][j]+(1.f-frac)*t[i][j-1]; 01047 n++; 01048 if(n == ly) 01049 break; 01050 ASSERT( yy[n] > yy[n-1] ); 01051 yl = yy[n]; 01052 } 01053 } 01054 01055 return 0; 01056 } 01057 STATIC void fhunt(realnum *xx, long int n, realnum x, long int *j) 01058 { 01059 /*lint -e731 boolean argument to equal / not equal */ 01060 long int jl, jm, jh, in; 01061 int up; 01062 01063 jl = *j; 01064 up = (xx[n-1] > xx[0]); 01065 if(jl < 0 || jl >= n) 01066 { 01067 jl = -1; 01068 jh = n; 01069 } 01070 else 01071 { 01072 in = 1; 01073 if((x >= xx[jl]) == up) 01074 { 01075 if(jl == n-1) 01076 { 01077 *j = jl; 01078 return; 01079 } 01080 jh = jl + 1; 01081 while ((x >= xx[jh]) == up) 01082 { 01083 jl = jh; 01084 in += in; 01085 jh += in; 01086 if(jh >= n) 01087 { 01088 jh = n; 01089 break; 01090 } 01091 } 01092 } 01093 else 01094 { 01095 if(jl == 0) 01096 { 01097 *j = -1; 01098 return; 01099 } 01100 jh = jl--; 01101 while ((x < xx[jl]) == up) 01102 { 01103 jh = jl; 01104 in += in; 01105 jl -= in; 01106 if(jl <= 0) 01107 { 01108 jl = 0; 01109 break; 01110 } 01111 } 01112 } 01113 } 01114 while (jh-jl != 1) 01115 { 01116 jm = (jh+jl)/2; 01117 if((x > xx[jm]) == up) 01118 jl = jm; 01119 else 01120 jh = jm; 01121 } 01122 *j = jl; 01123 return; 01124 } 01125 /*lint +e731 boolean argument to equal / not equal */