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 #include "cddefines.h" 00004 #include "phycon.h" 00005 #include "taulines.h" 00006 #include "mole.h" 00007 #include "atoms.h" 00008 #include "string.h" 00009 #include "thirdparty.h" 00010 #include "dense.h" 00011 #include "conv.h" 00012 #include "h2.h" 00013 #include "physconst.h" 00014 00015 void states_popfill( void); 00016 STATIC double LeidenCollRate(long, long, long, long,double); 00017 STATIC double Chianti_Upsilon(long, long, long, long,double); 00018 00019 #define DEBUGSTATE false 00020 00021 00022 /*Solving for the level populations*/ 00023 00024 void atmol_popsolve(void ) 00025 { 00026 realnum abund; 00027 DEBUG_ENTRY( "atmol_popsolve()" ); 00028 00029 for( long i=0; i<nSpecies; i++ ) 00030 { 00031 const char *spName = Species[i].chptrSpName; 00032 double *g, *ex, *pops, *depart; 00033 double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate,fcolldensity[9]; 00034 double *source, *sink; 00035 double cooltl, coolder; 00036 double fupsilon; 00037 long ipLo,ipHi,intCollNo; 00038 int nNegPop; 00039 bool lgZeroPop, lgDeBug = false; 00040 long nlev = Species[i].numLevels_max; 00041 00042 /* first find current density (cm-3) of species */ 00043 if( lgSpeciesMolecule[i] ) 00044 { 00045 struct molecule *SpeciesCurrent; 00048 if( (SpeciesCurrent = findspecies(Species[i].chptrSpName)) == &null_mole ) 00049 { 00050 /* did not find the species - print warning for now */ 00051 fprintf(ioQQQ," PROBLEM atmol_popsolve did not find molecular species %li\n",i); 00052 } 00053 abund = SpeciesCurrent->hevmol; 00054 } 00055 else 00056 { 00057 /* an atom or ion */ 00058 ASSERT( Species[i].intAtNo<LIMELM && Species[i].intIonStage<=Species[i].intAtNo); 00059 abund = dense.xIonDense[Species[i].intAtNo][Species[i].intIonStage]; 00060 } 00061 00062 if(abund<=SMALLFLOAT) 00063 { 00064 /* zero out populations and intensities */ 00065 /* NB - can't use states_popfill here, because it zeros ALL species */ 00066 atmolStates[i][0].Pop = 0.; 00067 for(long ipHi = 1; ipHi < nlev; ipHi++ ) 00068 { 00069 atmolStates[i][ipHi].Pop = 0.; 00070 for(long ipLo = 0; ipLo < ipHi; ipLo++ ) 00071 { 00072 if( atmolTrans[i][ipHi][ipLo].ipCont > 0 ) 00073 { 00074 atmolTrans[i][ipHi][ipLo].Emis->phots = 0.; 00075 atmolTrans[i][ipHi][ipLo].Emis->xIntensity = 0.; 00076 } 00077 } 00078 } 00079 continue; 00080 } 00081 00082 g = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00083 ex = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00084 pops = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00085 depart = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00086 source = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00087 sink = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00088 00089 AulEscp = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *)); 00090 col_str = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *)); 00091 AulDest = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *)); 00092 AulPump = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *)); 00093 CollRate = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *)); 00094 00095 for( long j=0; j< nlev; j++ ) 00096 { 00097 AulEscp[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00098 col_str[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00099 AulDest[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00100 AulPump[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00101 CollRate[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 00102 } 00103 00104 for( ipLo = 0; ipLo < nlev; ipLo++ ) 00105 { 00106 /*Filling in statistical weights & Excitation Energies*/ 00107 g[ipLo] = atmolStates[i][ipLo].g ; 00108 ex[ipLo] = atmolStates[i][ipLo].energy; 00109 /* zero some rates */ 00110 source[ipLo] = 0.; 00111 sink[ipLo] = 0.; 00112 AulEscp[ipLo][ipLo] = 0.; 00113 AulDest[ipLo][ipLo] = 0.; 00114 AulPump[ipLo][ipLo] = 0.; 00115 } 00116 00117 for( ipHi=1; ipHi<nlev; ipHi++ ) 00118 { 00119 for( ipLo=0; ipLo<ipHi; ipLo++ ) 00120 { 00121 if( atmolTrans[i][ipHi][ipLo].ipCont > 0 ) 00122 { 00123 AulEscp[ipHi][ipLo] = atmolTrans[i][ipHi][ipLo].Emis->Aul* 00124 (atmolTrans[i][ipHi][ipLo].Emis->Pesc + 00125 atmolTrans[i][ipHi][ipLo].Emis->Pelec_esc); 00126 AulDest[ipHi][ipLo] = atmolTrans[i][ipHi][ipLo].Emis->Aul* 00127 atmolTrans[i][ipHi][ipLo].Emis->Pdest; 00128 AulPump[ipLo][ipHi] = atmolTrans[i][ipHi][ipLo].Emis->pump; 00129 } 00130 else 00131 { 00132 AulEscp[ipHi][ipLo] = SMALLFLOAT; 00133 AulDest[ipHi][ipLo] = SMALLFLOAT; 00134 AulPump[ipLo][ipHi] = SMALLFLOAT; 00135 } 00136 00137 AulEscp[ipLo][ipHi] = 0.; 00138 AulDest[ipLo][ipHi] = 0.; 00139 AulPump[ipHi][ipLo] = 0.; 00140 } 00141 } 00142 00143 /*Setting all the collision strengths and collision rate to zero*/ 00144 for( ipHi= 0; ipHi<nlev; ipHi++) 00145 { 00146 for( ipLo= 0; ipLo<nlev; ipLo++) 00147 { 00148 col_str[ipHi][ipLo] = 0.; 00149 CollRate[ipHi][ipLo] = 0.; 00150 } 00151 } 00152 00153 /*Updating the CollRatesArray*/ 00154 /*You need to do this for a species indexed by i*/ 00155 for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++) 00156 { 00157 for( ipHi=1; ipHi<nlev; ipHi++) 00158 { 00159 for( ipLo=0; ipLo<ipHi; ipLo++) 00160 { 00161 /* molecule */ 00162 if( lgSpeciesMolecule[i] ) 00163 { 00164 if(CollRatesArray[i][intCollNo]!=NULL) 00165 { 00166 /*using the collision rate coefficients directly*/ 00167 CollRatesArray[i][intCollNo][ipHi][ipLo] = 00168 LeidenCollRate(i, intCollNo, ipHi, ipLo, phycon.te); 00169 } 00170 } 00171 /* atom or ion */ 00172 else 00173 { 00174 if(CollRatesArray[i][intCollNo]!=NULL) 00175 { 00176 fupsilon = Chianti_Upsilon(i, intCollNo, ipHi, ipLo, phycon.te); 00177 /*converting the collision strength to a collision rate coefficient*/ 00178 /*This formula converting collision strength to collision rate coefficent works fine for the electrons*/ 00179 /*For any other collider the mass would be different*/ 00180 CollRatesArray[i][intCollNo][ipHi][ipLo] = ((8.6291e-6)*fupsilon)/(g[ipHi]*phycon.sqrte); 00181 } 00182 } 00183 00184 /* now get the corresponding excitation rate */ 00185 if(CollRatesArray[i][intCollNo]!=NULL) 00186 { 00187 CollRatesArray[i][intCollNo][ipLo][ipHi] = 00188 CollRatesArray[i][intCollNo][ipHi][ipLo] * g[ipHi] / g[ipLo] * 00189 sexp( atmolTrans[i][ipHi][ipLo].EnergyK / phycon.te ); 00190 } 00191 } 00192 } 00193 } 00194 /*Get the densities seperately*/ 00195 /*Electron*/ 00196 fcolldensity[0] = dense.eden; 00197 /*Proton*/ 00198 fcolldensity[1] = dense.xIonDense[1][1]; 00199 /*Atomic Hydrogen*/ 00200 fcolldensity[2] = dense.xIonDense[1][0]; 00201 /*Helium*/ 00202 fcolldensity[3] = dense.xIonDense[2][0]; 00203 /*Helium+*/ 00204 fcolldensity[4] = dense.xIonDense[2][1]; 00205 /*Helium ++*/ 00206 fcolldensity[5] = dense.xIonDense[2][2]; 00207 /*Molecular Hydrogen*/ 00208 fcolldensity[8] = findspecies("H2")->hevmol; 00209 /*Ortho(6) and Para(7) Mol Hydrogen*/ 00210 if(h2.lgH2ON) 00211 { 00212 fcolldensity[6] = h2.ortho_density; 00213 fcolldensity[7] = h2.para_density; 00214 } 00215 else 00216 { 00217 fcolldensity[6] = (fcolldensity[8])*(0.75); 00218 fcolldensity[7] = (fcolldensity[8])*(0.25); 00219 } 00220 00221 /*Updating the CollRate*/ 00222 for( ipHi=0; ipHi<nlev; ipHi++ ) 00223 { 00224 for( ipLo=0; ipLo<nlev; ipLo++ ) 00225 { 00226 if( ipHi == ipLo ) 00227 { 00228 CollRate[ipHi][ipLo] = 0.; 00229 continue; 00230 } 00231 00232 for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++ ) 00233 { 00234 if(CollRatesArray[i][intCollNo]!=NULL) 00235 { 00236 CollRate[ipHi][ipLo] += CollRatesArray[i][intCollNo][ipHi][ipLo]*fcolldensity[intCollNo]; 00237 } 00238 } 00239 /*Correction for Helium*/ 00240 /*To include the effects of collision with Helium,the user must multiply the density by 1.14*/ 00241 /*pg 5,Schoier et al 2004*/ 00242 if(lgSpeciesMolecule[i]) 00243 { 00244 /*The collision rate coefficients for helium should not be present and that for molecular hydrogen should be present*/ 00245 if(CollRatesArray[i][3]==NULL && CollRatesArray[i][8]!=NULL ) 00246 { 00247 CollRate[ipHi][ipLo] += 0.7 *CollRatesArray[i][8][ipHi][ipLo]*fcolldensity[3]; 00248 } 00249 } 00250 } 00251 } 00253 /*Level Lowering as seen in mole_co_atom.cpp*/ 00254 /*Excitation energies need to be converted to temperature*/ 00255 /*ex[nUsed]/phycon.te cannot be much larger than 25, or matrix inversion will fail*/ 00256 while ( ((ex[nlev-1]*T1CM) > phycon.te*25.) && (nlev > 1) ) 00257 { 00258 --nlev; 00259 } 00260 Species[i].numLevels_local = nlev; 00261 00262 /* build some arrays */ 00263 atom_levelN( 00264 /* nlev is the number of levels to compute*/ 00265 nlev, 00266 /* ABUND is total abundance of species, used for nth equation 00267 * if balance equations are homogeneous */ 00268 abund, 00269 /* G(nlev) is stat weight of levels */ 00270 g, 00271 /* EX(nlev) is excitation potential of levels, deg K or wavenumbers 00272 * 0 for lowest level, all are energy rel to ground NOT d(ENER) */ 00273 ex, 00274 /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */ 00275 'w', 00276 /* populations [cm-3] of each level as deduced here */ 00277 pops, 00278 /* departure coefficient, derived below */ 00279 depart, 00280 /* net transition rate, A * esc prob, s-1 */ 00281 &AulEscp, 00282 /* col str from up to low */ 00283 &col_str, 00284 /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob, 00285 * asserts confirm that [ihi][ilo] is zero */ 00286 &AulDest, 00287 /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero */ 00288 &AulPump, 00289 /* collision rates (s^-1), evaluated here and returned for cooling by calling function, 00290 * unless following flag is true. If true then calling function has already filled 00291 * in these rates. CollRate[i][j] is rate from i to j */ 00292 &CollRate, 00293 /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */ 00294 source, 00295 /* this is an additional destruction rate to continuum, normally zero, units s-1 */ 00296 sink, 00297 /* flag saying whether CollRate already done (true), or we need to do it here (false), 00298 * this is stored in data)[ihi][ilo] as either downward rate or collis strength*/ 00299 true, 00300 /* total cooling and its derivative, set here but nothing done with it*/ 00301 &cooltl, 00302 &coolder, 00303 /* string used to identify calling program in case of error */ 00304 spName, 00305 /* nNegPop flag indicating what we have done 00306 * positive if negative populations occurred 00307 * zero if normal calculation done 00308 * negative if too cold (for some atoms other routine will be called in this case) */ 00309 &nNegPop, 00310 /* true if populations are zero, either due to zero abundance of very low temperature */ 00311 &lgZeroPop, 00312 /* option to print debug information */ 00313 lgDeBug ); 00314 00315 if( nNegPop == 0 ) 00316 { 00317 ASSERT( lgZeroPop == false ); 00318 00319 for( long j=0; j< nlev; j++ ) 00320 { 00321 atmolStates[i][j].Pop = pops[j]; 00322 } 00323 } 00324 else if( nNegPop > 0 ) 00325 { 00326 /* something should be done here. */ 00327 /* try another solver? */ 00328 /* probably don't want to quit */ 00329 fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population .\n"); 00330 cdEXIT( EXIT_FAILURE ); 00331 } 00332 else 00333 { 00334 /*In search phase we dont lower the levels*/ 00335 if(conv.lgSearch!= true) 00336 { 00337 for(long j=Species[i].numLevels_max-1; j>0; j--) 00338 { 00339 if((atmolStates[i][j].Pop/abund)<POPTHRES) 00340 { 00341 nlev--; 00342 } 00343 else 00344 break; 00345 } 00346 00347 ASSERT( nlev >= 1 ); 00348 /*Setting the number of energy levels to new limit*/ 00349 Species[i].numLevels_local = nlev; 00350 00351 if( nlev == 1 ) 00352 { 00353 atmolStates[i][0].Pop = abund; 00354 for( long j=1; j< Species[i].numLevels_max; j++ ) 00355 { 00356 atmolStates[i][j].Pop = 0.; 00357 } 00358 } 00359 else 00360 { 00361 /*Solving for level populations*/ 00362 atom_levelN(nlev,abund,g,ex,'w',pops,depart,&AulEscp,&col_str, 00363 &AulDest, &AulPump, &CollRate, source, sink,true,&cooltl, 00364 &coolder, spName, &nNegPop, &lgZeroPop, lgDeBug ); 00365 00366 for( long j=0; j< nlev; j++ ) 00367 { 00368 atmolStates[i][j].Pop = pops[j]; 00369 } 00370 for( long j=nlev; j< Species[i].numLevels_max; j++ ) 00371 { 00372 atmolStates[i][j].Pop = 0.; 00373 } 00374 } 00375 } 00376 } 00377 00378 /*Atmol line*/ 00379 for(long ipHi = 1; ipHi < nlev; ipHi++ ) 00380 { 00381 for(long ipLo = 0; ipLo < ipHi; ipLo++ ) 00382 { 00383 if( atmolTrans[i][ipHi][ipLo].ipCont > 0 ) 00384 { 00385 atmolTrans[i][ipHi][ipLo].Emis->phots = 00386 atmolTrans[i][ipHi][ipLo].Emis->Aul * 00387 atmolTrans[i][ipHi][ipLo].Hi->Pop * 00388 (atmolTrans[i][ipHi][ipLo].Emis->Pesc + 00389 atmolTrans[i][ipHi][ipLo].Emis->Pelec_esc); 00390 00391 atmolTrans[i][ipHi][ipLo].Emis->xIntensity = 00392 atmolTrans[i][ipHi][ipLo].Emis->phots * 00393 atmolTrans[i][ipHi][ipLo].EnergyErg; 00394 00395 /* population of lower level rel to ion, corrected for stim em */ 00396 00397 atmolTrans[i][ipHi][ipLo].Emis->PopOpc = 00398 atmolTrans[i][ipHi][ipLo].Lo->Pop - atmolTrans[i][ipHi][ipLo].Hi->Pop* 00399 atmolTrans[i][ipHi][ipLo].Lo->g/atmolTrans[i][ipHi][ipLo].Hi->g; 00400 fixit(); /* need to do something with these! */ 00401 atmolTrans[i][ipHi][ipLo].Emis->ColOvTot = 0.; 00402 atmolTrans[i][ipHi][ipLo].Coll.cool = 0.; 00403 atmolTrans[i][ipHi][ipLo].Coll.heat = 0.; 00404 } 00405 } 00406 } 00407 00408 /* option to print departure coefficients */ 00409 { 00410 enum {DEBUG_LOC=false}; 00411 00412 if( DEBUG_LOC ) 00413 { 00414 fprintf( ioQQQ, " Departure coefficients for species %li\n", i ); 00415 for( long j=0; j< nlev; j++ ) 00416 { 00417 fprintf( ioQQQ, " level %li \t Depar Coef %e\n", j, depart[j] ); 00418 } 00419 } 00420 } 00421 00422 /* free vectors */ 00423 free( g ); 00424 free( ex ); 00425 free( pops ); 00426 free( depart ); 00427 free( source ); 00428 free( sink ); 00429 00430 /* free arrays */ 00431 for( long j=0; j< nlev; j++ ) 00432 { 00433 free( AulEscp[j] ); 00434 free( col_str[j] ); 00435 free( AulDest[j] ); 00436 free( AulPump[j] ); 00437 free( CollRate[j] ); 00438 } 00439 free( AulEscp ); 00440 free( col_str ); 00441 free( AulDest ); 00442 free( AulPump ); 00443 free( CollRate ); 00444 } 00445 return; 00446 } 00447 00448 /*This function fills the population of the states to a dummy value of 0*/ 00449 void states_popfill( void) 00450 { 00451 DEBUG_ENTRY( "states_popfill()" ); 00452 00453 for( long i=0; i<nSpecies; i++) 00454 { 00455 for( long j=0; j<Species[i].numLevels_max; j++) 00456 { 00457 atmolStates[i][j].Pop = 0.; 00458 } 00459 } 00460 return; 00461 } 00462 00463 /*Leiden*/ 00464 double LeidenCollRate(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp) 00465 { 00466 double ret_collrate; 00467 int inttemps; 00468 DEBUG_ENTRY("LeidenCollRate()"); 00469 inttemps = AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps; 00470 00471 if( AtmolCollRateCoeff[ipSpecies][ipCollider].temps == NULL ) 00472 { 00473 return 0.; 00474 } 00475 00476 if(ftemp < AtmolCollRateCoeff[ipSpecies][ipCollider].temps[0]) 00477 { 00478 ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][0]; 00479 } 00480 else if(ftemp > AtmolCollRateCoeff[ipSpecies][ipCollider].temps[inttemps-1]) 00481 { 00482 ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][inttemps-1]; 00483 } 00484 else 00485 { 00486 ret_collrate = linint(AtmolCollRateCoeff[ipSpecies][ipCollider].temps, 00487 AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo], 00488 AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps, 00489 ftemp); 00490 } 00491 00492 if(DEBUGSTATE) 00493 printf("\nThe collision rate at temperature %f is %e",ftemp,ret_collrate); 00494 00495 ASSERT( !isnan( ret_collrate ) ); 00496 return(ret_collrate); 00497 00498 } 00499 00500 /*Chianti*/ 00501 double Chianti_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp) 00502 { 00503 double fdeltae,fscalingparam,fkte,fxt,fsups,fups; 00504 double *xs,*spl,*spl2; 00505 int i,intxs,inttype,intsplinepts; 00506 00507 DEBUG_ENTRY("Chianti_Upsilon()"); 00508 00509 if( AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL ) 00510 { 00511 return 0.; 00512 } 00513 00514 intsplinepts = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].nSplinePts; 00515 inttype = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].intTranType; 00516 fdeltae = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].EnergyDiff; 00517 fscalingparam = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].ScalingParam; 00518 00519 fkte = ftemp/fdeltae/1.57888e5; 00520 00521 /*Way the temperature is scaled*/ 00522 /*Burgess&Tully 1992:Paper gives only types 1 to 4*/ 00523 /*Found that the expressions were the same for 5 & 6 from the associated routine DESCALE_ALL*/ 00524 /*What about 7,8&9?*/ 00525 if(inttype ==1 || inttype == 4) 00526 { 00527 fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam))); 00528 } 00529 else if(inttype == 2 || inttype == 3||inttype == 5 || inttype == 6) 00530 { 00531 fxt = fkte/(fkte+fscalingparam); 00532 } 00533 else 00534 TotalInsanity(); 00535 00536 /*Creating spline points array*/ 00537 if(intsplinepts == 5) 00538 { 00539 xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00540 spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00541 spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00542 for(intxs=0;intxs<5;intxs++) 00543 { 00544 xs[intxs] = 0.25*intxs; 00545 spl[intxs] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline[intxs]; 00546 if(DEBUGSTATE) 00547 { 00548 printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]); 00549 getchar(); 00550 } 00551 } 00552 } 00553 else if(intsplinepts == 9) 00554 { 00555 xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00556 spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00557 spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00558 for( intxs=0; intxs<9; intxs++ ) 00559 { 00560 xs[intxs] = 0.125*intxs; 00561 spl[intxs] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline[intxs]; 00562 if(DEBUGSTATE) 00563 { 00564 printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]); 00565 getchar(); 00566 } 00567 } 00568 } 00569 else 00570 { 00571 TotalInsanity(); 00572 } 00573 00574 /*Finding the second derivative*/ 00575 for( i=0; i<intsplinepts; i++) 00576 { 00577 spl2[i] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].SplineSecDer[i]; 00578 } 00579 00580 if(DEBUGSTATE) 00581 { 00582 printf("\n"); 00583 for(intxs=0;intxs<intsplinepts;intxs++) 00584 { 00585 printf("The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]); 00586 } 00587 } 00588 00589 /*Extracting out the value*/ 00590 splint(xs,spl,spl2,intsplinepts,fxt,&fsups); 00591 00592 /*Finding upsilon*/ 00593 if(inttype == 1) 00594 { 00595 fups = fsups*log(fkte+exp(1.0)); 00596 } 00597 else if(inttype == 2) 00598 { 00599 fups = fsups; 00600 } 00601 else if(inttype == 3) 00602 { 00603 fups = fsups/(fkte+1.0) ; 00604 } 00605 else if(inttype == 4) 00606 { 00607 fups = fsups*log(fkte+fscalingparam) ; 00608 } 00609 else if(inttype == 5) 00610 { 00611 fups = fsups/fkte ; 00612 } 00613 else if(inttype == 6) 00614 { 00615 fups = pow(10.0,fsups) ; 00616 } 00617 else 00618 { 00619 TotalInsanity(); 00620 } 00621 00622 free( xs ); 00623 free( spl ); 00624 free( spl2 ); 00625 00626 ASSERT(fups>=0); 00627 return(fups); 00628 }