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