cloudy trunk

temp_change.cpp

Go to the documentation of this file.
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 */
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated for cloudy by doxygen 1.7.3