00001
00002
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 #include "secondaries.h"
00015 #include "thermal.h"
00016 #include "cooling.h"
00017 #include "lines_service.h"
00018
00019 void states_popfill( void);
00020 STATIC double LeidenCollRate(long, long, long, long,double);
00021 STATIC double CHIANTI_Upsilon(long, long, long, long,double);
00022
00023 #define DEBUGSTATE false
00024
00025 static double *g, *ex, *pops, *depart, *source, *sink;
00026 static double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate;
00027
00028
00029
00030 void dBase_solve(void )
00031 {
00032 realnum abund;
00033 DEBUG_ENTRY( "dBase_solve()" );
00034 static bool lgFirstPass = true;
00035 static long maxNumLevels = 1;
00036 double totalHeating = 0.;
00037
00038 if( nSpecies==0 )
00039 return;
00040
00041 if( lgFirstPass )
00042 {
00043 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
00044 maxNumLevels = MAX2( maxNumLevels, Species[ipSpecies].numLevels_max );
00045
00046 g = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00047 ex = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00048 pops = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00049 depart = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00050 source = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00051 sink = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00052
00053 AulEscp = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
00054 col_str = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
00055 AulDest = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
00056 AulPump = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
00057 CollRate = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
00058
00059 for( long j=0; j< maxNumLevels; j++ )
00060 {
00061 AulEscp[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00062 col_str[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00063 AulDest[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00064 AulPump[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00065 CollRate[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
00066 }
00067
00068 lgFirstPass = false;
00069 }
00070
00071
00072 memset( g, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00073 memset( ex, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00074 memset( pops, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00075 memset( depart, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00076 memset( source, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00077 memset( sink, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00078 for( long j=0; j< maxNumLevels; j++ )
00079 {
00080 memset( AulEscp[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00081 memset( col_str[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00082 memset( AulDest[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00083 memset( AulPump[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00084 memset( CollRate[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
00085 }
00086
00087
00088 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
00089 {
00090 const char *spName = Species[ipSpecies].chLabel;
00091 double fcolldensity[9];
00092 double cooltl, coolder;
00093 double fupsilon;
00094 long intCollNo;
00095 int nNegPop;
00096 bool lgZeroPop, lgDeBug = false;
00097 double totaldCdT = 0.;
00098
00099 Species[ipSpecies].CoolTotal = 0.;
00100
00101 #if 0
00102
00103 Species[ipSpecies].numLevels_local = MIN2( Species[ipSpecies].numLevels_local, 10 );
00104 #endif
00105
00106
00107 if( Species[ipSpecies].lgMolecular )
00108 {
00109 struct molecule *SpeciesCurrent;
00112 if( (SpeciesCurrent = findspecies(Species[ipSpecies].chLabel)) == &null_mole )
00113 {
00114
00115 if( !conv.nTotalIoniz )
00116 fprintf(ioQQQ," PROBLEM dBase_solve did not find molecular species %li\n",ipSpecies);
00117 }
00118 abund = SpeciesCurrent->hevmol;
00119 }
00120 else
00121 {
00122
00123 ASSERT( dBaseStates[ipSpecies][0].nelem<=LIMELM && dBaseStates[ipSpecies][0].IonStg<=dBaseStates[ipSpecies][0].nelem+1 );
00124 abund = dense.xIonDense[ dBaseStates[ipSpecies][0].nelem-1 ][ dBaseStates[ipSpecies][0].IonStg-1 ];
00125 }
00126
00127 abund *= Species[ipSpecies].fracType * Species[ipSpecies].fracIsotopologue;
00128
00129
00130 if( conv.nTotalIoniz == 0)
00131 Species[ipSpecies].lgActive = true;
00132
00133 bool lgMakeInActive = (abund <= 1e-20 * dense.xNucleiTotal);
00134 if( lgMakeInActive && Species[ipSpecies].lgActive )
00135 {
00136
00137 dBaseStates[ipSpecies][0].Pop = 0.;
00138 for(long ipHi = 1; ipHi < Species[ipSpecies].numLevels_max; ipHi++ )
00139 {
00140 dBaseStates[ipSpecies][ipHi].Pop = 0.;
00141 transition *t = &dBaseTrans[ipSpecies][ipHi][0];
00142 for(long ipLo = 0; ipLo < ipHi; ipLo++ )
00143 {
00144 t->Emis->phots = 0.;
00145 t->Emis->xIntensity = 0.;
00146 t->Coll.col_str = 0.;
00147 t->Coll.cool = 0.;
00148 t->Coll.heat = 0.;
00149 ++t;
00150 }
00151 }
00152 Species[ipSpecies].lgActive = false;
00153 }
00154
00155 if( !lgMakeInActive )
00156 Species[ipSpecies].lgActive = true;
00157
00158 if( !Species[ipSpecies].lgActive )
00159 continue;
00160
00161
00162 if( conv.lgSearch )
00163 Species[ipSpecies].numLevels_local = Species[ipSpecies].numLevels_max;
00164
00165 for( long ipLo = 0; ipLo < Species[ipSpecies].numLevels_local; ipLo++ )
00166 {
00167
00168 g[ipLo] = dBaseStates[ipSpecies][ipLo].g ;
00169
00170
00171 ex[ipLo] = dBaseStates[ipSpecies][ipLo].energy.WN() -
00172 dBaseStates[ipSpecies][0].energy.WN();
00173
00174 source[ipLo] = 0.;
00175 sink[ipLo] = 0.;
00176 AulEscp[ipLo][ipLo] = 0.;
00177 AulDest[ipLo][ipLo] = 0.;
00178 AulPump[ipLo][ipLo] = 0.;
00179 }
00180
00181
00182 if( ex[0] <= dBaseStates[ipSpecies][0].energy.WN()* 10. *DBL_EPSILON )
00183 ex[0] = 0.;
00184 else
00185 TotalInsanity();
00186
00187 for( long ipHi=1; ipHi<Species[ipSpecies].numLevels_local; ipHi++ )
00188 {
00189 for( long ipLo=0; ipLo<ipHi; ipLo++ )
00190 {
00191 if( dBaseTrans[ipSpecies][ipHi][ipLo].ipCont > 0 )
00192 {
00193 AulEscp[ipHi][ipLo] = dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Aul*
00194 (dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Pesc +
00195 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Pelec_esc);
00196 AulDest[ipHi][ipLo] = dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Aul*
00197 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Pdest;
00198 AulPump[ipLo][ipHi] = dBaseTrans[ipSpecies][ipHi][ipLo].Emis->pump;
00199 }
00200 else
00201 {
00202 AulEscp[ipHi][ipLo] = SMALLFLOAT;
00203 AulDest[ipHi][ipLo] = SMALLFLOAT;
00204 AulPump[ipLo][ipHi] = SMALLFLOAT;
00205 }
00206
00207 AulEscp[ipLo][ipHi] = 0.;
00208 AulDest[ipLo][ipHi] = 0.;
00209 AulPump[ipHi][ipLo] = 0.;
00210 }
00211 }
00212
00213
00214 for( long ipHi= 0; ipHi<Species[ipSpecies].numLevels_local; ipHi++)
00215 {
00216 for( long ipLo= 0; ipLo<Species[ipSpecies].numLevels_local; ipLo++)
00217 {
00218 col_str[ipHi][ipLo] = 0.;
00219 CollRate[ipHi][ipLo] = 0.;
00220 }
00221 }
00222
00223
00224
00225 for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++)
00226 {
00227 for( long ipHi=1; ipHi<Species[ipSpecies].numLevels_local; ipHi++)
00228 {
00229 for( long ipLo=0; ipLo<ipHi; ipLo++)
00230 {
00231
00232 if( Species[ipSpecies].lgMolecular )
00233 {
00234 if(CollRatesArray[ipSpecies][intCollNo]!=NULL)
00235 {
00236
00237 CollRatesArray[ipSpecies][intCollNo][ipHi][ipLo] =
00238 LeidenCollRate(ipSpecies, intCollNo, ipHi, ipLo, phycon.te);
00239 }
00240 }
00241
00242 else
00243 {
00244 if(CollRatesArray[ipSpecies][intCollNo]!=NULL)
00245 {
00246 fupsilon = CHIANTI_Upsilon(ipSpecies, intCollNo, ipHi, ipLo, phycon.te);
00247
00248
00249
00250 if( intCollNo==1 && AtmolCollSplines[ipSpecies][ipHi][ipLo][intCollNo].intTranType==6 )
00251 {
00252 CollRatesArray[ipSpecies][intCollNo][ipHi][ipLo] = fupsilon;
00253 }
00254 else
00255 {
00256
00257
00258
00259 ASSERT( intCollNo == 0 );
00260 CollRatesArray[ipSpecies][intCollNo][ipHi][ipLo] = (COLL_CONST*fupsilon)/(g[ipHi]*phycon.sqrte);
00261 }
00262 }
00263 }
00264
00265
00266 if(CollRatesArray[ipSpecies][intCollNo]!=NULL)
00267 {
00268 CollRatesArray[ipSpecies][intCollNo][ipLo][ipHi] =
00269 CollRatesArray[ipSpecies][intCollNo][ipHi][ipLo] * g[ipHi] / g[ipLo] *
00270 sexp( dBaseTrans[ipSpecies][ipHi][ipLo].EnergyK / phycon.te );
00271 }
00272 }
00273 }
00274 }
00275
00276
00277 fcolldensity[0] = dense.eden;
00278
00279 fcolldensity[1] = dense.xIonDense[ipHYDROGEN][1];
00280
00281 fcolldensity[2] = dense.xIonDense[ipHYDROGEN][0];
00282
00283 fcolldensity[3] = dense.xIonDense[ipHELIUM][0];
00284
00285 fcolldensity[4] = dense.xIonDense[ipHELIUM][1];
00286
00287 fcolldensity[5] = dense.xIonDense[ipHELIUM][2];
00288
00289 fcolldensity[8] = findspecies("H2")->hevmol;
00290
00291 if(h2.lgH2ON)
00292 {
00293 fcolldensity[6] = h2.ortho_density;
00294 fcolldensity[7] = h2.para_density;
00295 }
00296 else
00297 {
00298 fcolldensity[6] = (fcolldensity[8])*(0.75);
00299 fcolldensity[7] = (fcolldensity[8])*(0.25);
00300 }
00301
00302
00303 for( long ipHi=0; ipHi<Species[ipSpecies].numLevels_local; ipHi++ )
00304 {
00305 for( long ipLo=0; ipLo<Species[ipSpecies].numLevels_local; ipLo++ )
00306 {
00307 if( ipHi == ipLo )
00308 {
00309 CollRate[ipHi][ipLo] = 0.;
00310 continue;
00311 }
00312
00313 for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++ )
00314 {
00315 if(CollRatesArray[ipSpecies][intCollNo]!=NULL)
00316 {
00317 CollRate[ipHi][ipLo] += CollRatesArray[ipSpecies][intCollNo][ipHi][ipLo]*fcolldensity[intCollNo];
00318 }
00319 }
00320
00321
00322
00323 if( Species[ipSpecies].lgMolecular )
00324 {
00325
00326 if(CollRatesArray[ipSpecies][3]==NULL && CollRatesArray[ipSpecies][8]!=NULL )
00327 {
00328 CollRate[ipHi][ipLo] += 0.7 *CollRatesArray[ipSpecies][8][ipHi][ipLo]*fcolldensity[3];
00329 }
00330
00331
00332 if(CollRatesArray[ipSpecies][2]==NULL )
00333 {
00334 double colliderDensity = fcolldensity[2];
00335 if( CollRatesArray[ipSpecies][3]!=NULL )
00336 {
00337 CollRate[ipHi][ipLo] += 2.0 *CollRatesArray[ipSpecies][3][ipHi][ipLo]*colliderDensity;
00338 }
00339 else if( CollRatesArray[ipSpecies][6]!=NULL )
00340 {
00341 CollRate[ipHi][ipLo] += 1.4 *CollRatesArray[ipSpecies][6][ipHi][ipLo]*colliderDensity;
00342 }
00343 else if( CollRatesArray[ipSpecies][7]!=NULL )
00344 {
00345 CollRate[ipHi][ipLo] += 1.4 *CollRatesArray[ipSpecies][7][ipHi][ipLo]*colliderDensity;
00346 }
00347 else if( CollRatesArray[ipSpecies][8]!=NULL )
00348 {
00349 CollRate[ipHi][ipLo] += 1.4 *CollRatesArray[ipSpecies][8][ipHi][ipLo]*colliderDensity;
00350 }
00351 else
00352 {
00353
00354 CollRate[ipHi][ipLo] += 1e-13 * colliderDensity;
00355 }
00356 }
00357
00358
00359 if(CollRatesArray[ipSpecies][1]==NULL )
00360 {
00361 double colliderDensity = fcolldensity[1];
00362 if( CollRatesArray[ipSpecies][4]!=NULL )
00363 {
00364 CollRate[ipHi][ipLo] += 2.0 *CollRatesArray[ipSpecies][3][ipHi][ipLo]*colliderDensity;
00365 }
00366 else
00367 {
00368
00369 CollRate[ipHi][ipLo] += 1e-13 * colliderDensity;
00370 }
00371 }
00372 }
00373
00374
00375 if( ipHi>ipLo )
00376 {
00377 fixit();
00378
00379
00380
00381
00382 if( dBaseTrans[ipSpecies][ipHi][ipLo].ipCont > 0 )
00383 {
00384
00385
00386
00387 double RateUp = secondaries.x12tot
00388
00389
00390
00391 ;
00392
00393 double RateDown = RateUp * dBaseTrans[ipSpecies][ipHi][ipLo].Lo->g /
00394 dBaseTrans[ipSpecies][ipHi][ipLo].Hi->g;
00395
00396 CollRate[ipLo][ipHi] += RateUp;
00397 CollRate[ipHi][ipLo] += RateDown;
00398 }
00399 }
00400 }
00401 }
00402
00403
00404 atom_levelN(
00405
00406 Species[ipSpecies].numLevels_local,
00407
00408
00409 abund,
00410
00411 g,
00412
00413
00414 ex,
00415
00416 'w',
00417
00418 pops,
00419
00420 depart,
00421
00422 &AulEscp,
00423
00424 &col_str,
00425
00426
00427 &AulDest,
00428
00429 &AulPump,
00430
00431
00432
00433 &CollRate,
00434
00435 source,
00436
00437 sink,
00438
00439
00440 true,
00441
00442 &cooltl,
00443 &coolder,
00444
00445 spName,
00446
00447
00448
00449
00450 &nNegPop,
00451
00452 &lgZeroPop,
00453
00454 lgDeBug );
00455
00456 if( nNegPop > 0 )
00457 {
00458
00459 fprintf(ioQQQ," PROBLEM in dBase_solve, atom_levelN returned negative population .\n");
00460 cdEXIT( EXIT_FAILURE );
00461 }
00462
00463
00464 while( (pops[Species[ipSpecies].numLevels_local-1]<=0 ) &&
00465 (Species[ipSpecies].numLevels_local > 1) )
00466 --Species[ipSpecies].numLevels_local;
00467
00468 for( long j=0;j< Species[ipSpecies].numLevels_local; j++ )
00469 dBaseStates[ipSpecies][j].Pop = MAX2(pops[j],SMALLFLOAT);
00470
00471 for( long j=Species[ipSpecies].numLevels_local;
00472 j< Species[ipSpecies].numLevels_max; j++ )
00473 dBaseStates[ipSpecies][j].Pop = 0.;
00474
00475
00476 for(long ipHi = 1; ipHi < Species[ipSpecies].numLevels_local; ipHi++ )
00477 {
00478 for(long ipLo = 0; ipLo < ipHi; ipLo++ )
00479 {
00480 double qLU = dBaseTrans[ipSpecies][ipHi][ipLo].Lo->Pop * CollRate[ipLo][ipHi];
00481 double qUL = dBaseTrans[ipSpecies][ipHi][ipLo].Hi->Pop * CollRate[ipHi][ipLo];
00482 double netCooling = (qLU - qUL) * dBaseTrans[ipSpecies][ipHi][ipLo].EnergyErg;
00483
00484 if( netCooling > 0. )
00485 {
00486 dBaseTrans[ipSpecies][ipHi][ipLo].Coll.cool = netCooling;
00487 dBaseTrans[ipSpecies][ipHi][ipLo].Coll.heat = 0.;
00488 }
00489 else
00490 {
00491 dBaseTrans[ipSpecies][ipHi][ipLo].Coll.cool = 0;
00492 dBaseTrans[ipSpecies][ipHi][ipLo].Coll.heat = -netCooling;
00493 }
00494
00495 Species[ipSpecies].CoolTotal += netCooling;
00496 totaldCdT += netCooling *
00497 (dBaseTrans[ipSpecies][ipHi][ipLo].EnergyK * thermal.tsq1 - thermal.halfte);
00498
00499 if( dBaseTrans[ipSpecies][ipHi][ipLo].ipCont > 0 )
00500 {
00501 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->phots =
00502 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Aul *
00503 dBaseTrans[ipSpecies][ipHi][ipLo].Hi->Pop *
00504 (dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Pesc +
00505 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->Pelec_esc);
00506
00507 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->xIntensity =
00508 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->phots *
00509 dBaseTrans[ipSpecies][ipHi][ipLo].EnergyErg;
00510
00511
00512 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->PopOpc =
00513 dBaseTrans[ipSpecies][ipHi][ipLo].Lo->Pop - dBaseTrans[ipSpecies][ipHi][ipLo].Hi->Pop*
00514 dBaseTrans[ipSpecies][ipHi][ipLo].Lo->g/dBaseTrans[ipSpecies][ipHi][ipLo].Hi->g;
00515
00516
00517 if( CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi] > 0. )
00518 {
00519 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->ColOvTot = CollRate[ipLo][ipHi]/
00520 (CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi]);
00521 }
00522 else
00523 dBaseTrans[ipSpecies][ipHi][ipLo].Emis->ColOvTot = 0.;
00524
00525
00526
00527
00528 dBaseTrans[ipSpecies][ipHi][ipLo].Coll.col_str = (realnum)(CollRate[ipHi][ipLo] *
00529 (g[ipHi]*phycon.sqrte)/COLL_CONST);
00530 }
00531 else
00532 dBaseTrans[ipSpecies][ipHi][ipLo].Coll.col_str = 0.;
00533 }
00534 }
00535
00536 fixit();
00537
00538
00539 for(long ipHi = Species[ipSpecies].numLevels_local; ipHi< Species[ipSpecies].numLevels_max; ipHi++ )
00540 {
00541 for(long ipLo = 0; ipLo < ipHi; ipLo++ )
00542 {
00543 EmLineZero( dBaseTrans[ipSpecies][ipHi][ipLo].Emis );
00544 }
00545 }
00546
00547 thermal.dCooldT += totaldCdT;
00548 CoolAdd( Species[ipSpecies].chLabel, 0. , MAX2(0.,Species[ipSpecies].CoolTotal) );
00549 totalHeating += MAX2(0., -Species[ipSpecies].CoolTotal);
00550
00551
00552 {
00553 enum {DEBUG_LOC=false};
00554
00555 if( DEBUG_LOC )
00556 {
00557 fprintf( ioQQQ, " Departure coefficients for species %li\n", ipSpecies );
00558 for( long j=0; j< Species[ipSpecies].numLevels_local; j++ )
00559 {
00560 fprintf( ioQQQ, " level %li \t Depar Coef %e\n", j, depart[j] );
00561 }
00562 }
00563 }
00564 }
00565
00566
00567 thermal.heating[0][27] = totalHeating;
00568
00569 return;
00570 }
00571
00572
00573 void states_popfill( void)
00574 {
00575 DEBUG_ENTRY( "states_popfill()" );
00576
00577 for( long i=0; i<nSpecies; i++)
00578 {
00579 for( long j=0; j<Species[i].numLevels_max; j++)
00580 {
00581 dBaseStates[i][j].Pop = 0.;
00582 }
00583 }
00584 return;
00585 }
00586
00587
00588 STATIC double LeidenCollRate(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
00589 {
00590 double ret_collrate;
00591 int inttemps;
00592 DEBUG_ENTRY("LeidenCollRate()");
00593 inttemps = AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps;
00594
00595 if( AtmolCollRateCoeff[ipSpecies][ipCollider].temps == NULL )
00596 {
00597 return 0.;
00598 }
00599
00600 if(ftemp < AtmolCollRateCoeff[ipSpecies][ipCollider].temps[0])
00601 {
00602 ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][0];
00603 }
00604 else if(ftemp > AtmolCollRateCoeff[ipSpecies][ipCollider].temps[inttemps-1])
00605 {
00606 ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][inttemps-1];
00607 }
00608 else if( AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps==1 )
00609 {
00610
00611 ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][0];
00612 }
00613 else
00614 {
00615 ret_collrate = linint(AtmolCollRateCoeff[ipSpecies][ipCollider].temps,
00616 AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo],
00617 AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps,
00618 ftemp);
00619 }
00620
00621 if(DEBUGSTATE)
00622 printf("\nThe collision rate at temperature %f is %e",ftemp,ret_collrate);
00623
00624 ASSERT( !isnan( ret_collrate ) );
00625 return(ret_collrate);
00626
00627 }
00628
00629
00630 STATIC double CHIANTI_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
00631 {
00632 double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
00633 int intxs,inttype,intsplinepts;
00634
00635 DEBUG_ENTRY("CHIANTI_Upsilon()");
00636
00637 if( AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL )
00638 {
00639 return 0.;
00640 }
00641
00642 intsplinepts = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].nSplinePts;
00643 inttype = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].intTranType;
00644 fdeltae = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].EnergyDiff;
00645 fscalingparam = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].ScalingParam;
00646
00647 fkte = ftemp/fdeltae/1.57888e5;
00648
00649
00650
00651
00652
00653 if( inttype ==1 || inttype==4 )
00654 {
00655 fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam)));
00656 }
00657 else if(inttype == 2 || inttype == 3||inttype == 5 || inttype == 6)
00658 {
00659 fxt = fkte/(fkte+fscalingparam);
00660 }
00661 else
00662 TotalInsanity();
00663
00664 double xs[9],*spl,*spl2;
00665
00666 if(intsplinepts == 5)
00667 {
00668 spl = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline;
00669 for(intxs=0;intxs<5;intxs++)
00670 {
00671 xs[intxs] = 0.25*intxs;
00672 if(DEBUGSTATE)
00673 {
00674 printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
00675 getchar();
00676 }
00677 }
00678 }
00679 else if(intsplinepts == 9)
00680 {
00681 spl = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline;
00682 for( intxs=0; intxs<9; intxs++ )
00683 {
00684 xs[intxs] = 0.125*intxs;
00685 if(DEBUGSTATE)
00686 {
00687 printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
00688 getchar();
00689 }
00690 }
00691 }
00692 else
00693 {
00694 TotalInsanity();
00695 }
00696
00697
00698 spl2 = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].SplineSecDer;
00699
00700 if(DEBUGSTATE)
00701 {
00702 printf("\n");
00703 for(intxs=0;intxs<intsplinepts;intxs++)
00704 {
00705 printf("The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]);
00706 }
00707 }
00708
00709
00710
00711
00712 fsups = linint( xs, spl, intsplinepts, fxt);
00713
00714
00715 if(inttype == 1)
00716 {
00717 fups = fsups*log(fkte+exp(1.0));
00718 }
00719 else if(inttype == 2)
00720 {
00721 fups = fsups;
00722 }
00723 else if(inttype == 3)
00724 {
00725 fups = fsups/(fkte+1.0) ;
00726 }
00727 else if(inttype == 4)
00728 {
00729 fups = fsups*log(fkte+fscalingparam) ;
00730 }
00731 else if(inttype == 5)
00732 {
00733 fups = fsups/fkte ;
00734 }
00735 else if(inttype == 6)
00736 {
00737 fups = pow(10.0,fsups) ;
00738 }
00739 else
00740 {
00741 TotalInsanity();
00742 }
00743
00744 if( fups < 0. )
00745 {
00746 fprintf( ioQQQ," WARNING: Negative upsilon in species %s, collider %li, indices %4li %4li, Te = %e.\n",
00747 Species[ipSpecies].chLabel, ipCollider, ipHi, ipLo, ftemp );
00748 fups = 0.;
00749 }
00750 ASSERT(fups>=0);
00751 return(fups);
00752 }
00753