cloudy trunk
Loading...
Searching...
No Matches
helike_cs.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*HeCollid evaluate collisional rates */
4/*HeCSInterp interpolate on He1 collision strengths */
5/*AtomCSInterp do the atom */
6/*IonCSInterp do the ions */
7/*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
8#include "cddefines.h"
9#include "atmdat.h"
10#include "conv.h"
11#include "dense.h"
12#include "helike.h"
13#include "helike_cs.h"
14#include "hydro_vs_rates.h"
15#include "iso.h"
16#include "lines_service.h"
17#include "opacity.h"
18#include "phycon.h"
19#include "physconst.h"
20#include "rfield.h"
21#include "taulines.h"
22#include "thirdparty.h"
23#include "trace.h"
24
26vector<double> CSTemp;
29
30/* returns thermally-averaged Seaton 62 collision strength. */
31STATIC double S62_Therm_ave_coll_str( double proj_energy_overKT, long nelem, long Collider, double deltaE, double osc_strength, double temp,
32 double stat_weight, double I_energy_eV );
33
34/* all of these are used in the calculation of Stark collision strengths
35 * following the algorithms in Vrinceanu & Flannery (2001). */
36STATIC double collision_strength_VF01( long ipISO, long nelem, long n, long l, long lp, long s, long Collider,
37 double ColliderCharge, double temp, double velOrEner, bool lgParamIsRedVel );
38STATIC double L_mix_integrand_VF01( long n, long l, long lp, double bmax, double red_vel, double an, double ColliderCharge, double alpha );
39STATIC double StarkCollTransProb_VF01( long int n, long int l, long int lp, double alpha, double deltaPhi);
40
41
43{
44public:
47
48 double operator() (double proj_energy_overKT)
49 {
50 double col_str = S62_Therm_ave_coll_str( proj_energy_overKT, nelem, Collider, deltaE, osc_strength,
52 return col_str;
53 }
54};
55
57{
58public:
59 long ipISO, nelem, n, l, lp, s, Collider;
62
63 double operator() (double EOverKT)
64 {
67 return exp( -1.*EOverKT ) * col_str;
68 }
69};
70
72{
73public:
74 long n, l, lp;
76
77 double operator() (double alpha)
78 {
79 double integrand = L_mix_integrand_VF01( n, l, lp,
80 bmax, red_vel, an, ColliderCharge, alpha );
81 return integrand;
82 }
83};
84
85
86/* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */
87static const double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
88static const double ColliderCharge[4] = {1.0, 1.0, 1.0, 2.0};
89
90void HeCollidSetup( void )
91{
92 /* this must be longer than data path string, set in path.h/cpu.cpp */
93 long i, i1, j, nelem, ipHi, ipLo;
94 bool lgEOL, lgHIT;
95 FILE *ioDATA;
96
97# define chLine_LENGTH 1000
98 char chLine[chLine_LENGTH];
99
100 DEBUG_ENTRY( "HeCollidSetup()" );
101
102 /* get the collision strength data for the He 1 lines */
103 if( trace.lgTrace )
104 fprintf( ioQQQ," HeCollidSetup opening he1_cs.dat:");
105
106 ioDATA = open_data( "he1_cs.dat", "r" );
107
108 /* check that magic number is ok */
109 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
110 {
111 fprintf( ioQQQ, " HeCollidSetup could not read first line of he1_cs.dat.\n");
113 }
114 i = 1;
115 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
116 /*i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
117 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);*/
118
119 /* the following is to check the numbers that appear at the start of he1_cs.dat */
120 if( i1 !=COLLISMAGIC )
121 {
122 fprintf( ioQQQ,
123 " HeCollidSetup: the version of he1_cs.dat is not the current version.\n" );
124 fprintf( ioQQQ,
125 " HeCollidSetup: I expected to find the number %i and got %li instead.\n" ,
126 COLLISMAGIC, i1 );
127 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
129 }
130
131 /* get the array of temperatures */
132 lgHIT = false;
133 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
134 {
135 /* only look at lines without '#' in first col */
136 if( chLine[0] == '#')
137 continue;
138
139 lgHIT = true;
140 lgEOL = false;
141 char *chTemp = strtok(chLine," \t\n");
142 while( chTemp != NULL )
143 {
144 CSTemp.push_back( atof(chTemp) );
145 chTemp = strtok(NULL," \t\n");
146 }
147 break;
148 }
149 if( !lgHIT )
150 {
151 fprintf( ioQQQ, " HeCollidSetup could not find line in CS temperatures in he1_cs.dat.\n");
153 }
154 ASSERT( CSTemp.size() == 9U );
155
156 /* create space for array of CS values, if not already done */
157 {
158 long nelem = ipHELIUM;
159 long numLevs = iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max;
160 HeCS.reserve( numLevs );
161 for( long ipHi=1; ipHi < numLevs; ++ipHi )
162 {
163 HeCS.reserve( ipHi, ipHi );
164 for( long ipLo=0; ipLo<ipHi; ++ipLo )
165 HeCS.reserve( ipHi, ipLo, CSTemp.size() );
166 }
167 HeCS.alloc();
168 }
169
170 /* now read in the CS data */
171 lgHIT = false;
172 nelem = ipHELIUM;
173 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
174 {
175 char *chTemp;
176 /* only look at lines without '#' in first col */
177 if( (chLine[0] == ' ') || (chLine[0]=='\n') )
178 break;
179 if( chLine[0] != '#')
180 {
181 lgHIT = true;
182
183 /* get lower and upper level index */
184 j = 1;
185 ipLo = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
186 ipHi = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
187 ASSERT( ipLo < ipHi );
188 if( ipHi >= iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
189 continue;
190 else
191 {
192 chTemp = chLine;
193 /* skip over 4 tabs to start of cs data */
194 for( long i=0; i<3; ++i )
195 {
196 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
197 {
198 fprintf( ioQQQ, " HeCollidSetup could not init cs\n" );
200 }
201 ++chTemp;
202 }
203
204 for( unsigned i=0; i< CSTemp.size(); ++i )
205 {
206 double a;
207 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
208 {
209 fprintf( ioQQQ, " HeCollidSetup could not scan cs, current indices: %li %li\n", ipHi, ipLo );
211 }
212 ++chTemp;
213 sscanf( chTemp , "%le" , &a );
214 HeCS[ipHi][ipLo][i] = (realnum)a;
215 }
216 }
217 }
218 }
219
220 /* close the data file */
221 fclose( ioDATA );
222
223 return;
224}
225
226/* Choose either AtomCSInterp or IonCSInterp */
227realnum HeCSInterp(long int nelem,
228 long int ipHi,
229 long int ipLo,
230 long int Collider )
231{
232 realnum cs, factor1;
233
234 /* This variable is for diagnostic purposes:
235 * a string used in the output to designate where each cs comes from. */
236 const char *where = " ";
237
238 DEBUG_ENTRY( "HeCSInterp()" );
239
240 if( !iso_ctrl.lgColl_excite[ipHE_LIKE] )
241 {
242 return (realnum)1E-10;
243 }
244
245 if( nelem == ipHELIUM )
246 {
247 /* do for helium */
248 cs = AtomCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
249 }
250 else
251 {
252 /* get collision strengths for an ion */
253 cs = IonCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
254 }
255
256 ASSERT( cs >= 0.f );
257
258 /* in many cases the correction factor for split states has already been made,
259 * if not then factor is still negative */
260 /* Remove the second test here when IonCSInterp is up to par with AtomCSInterp */
261 ASSERT( factor1 >= 0.f || nelem!=ipHELIUM );
262 if( factor1 < 0.f )
263 {
264 ASSERT( iso_ctrl.lgCS_Vriens[ipHE_LIKE] );
265
266 factor1 = 1.f;
267 }
268
269 /* take factor into account, usually 1, ratio of stat weights if within 2 3P
270 * and with collisions from collapsed to resolved levels */
271 cs *= factor1;
272
273 {
274 /*@-redef@*/
275 enum {DEBUG_LOC=false};
276 /*@+redef@*/
277
278 if( DEBUG_LOC && ( nelem==ipOXYGEN ) && (cs > 0.f) && (iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == 2)
279 && ( iso_sp[ipHE_LIKE][nelem].st[ipLo].n() <= 2 ) )
280 fprintf(ioQQQ,"%li\t%li\t%li\t-\t%li\t%li\t%li\t%.2e\t%s\n",
281 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(), iso_sp[ipHE_LIKE][nelem].st[ipLo].S() ,
282 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(), iso_sp[ipHE_LIKE][nelem].st[ipHi].n() ,
283 iso_sp[ipHE_LIKE][nelem].st[ipHi].S(), iso_sp[ipHE_LIKE][nelem].st[ipHi].l() , cs,where);
284 }
285
286 return MAX2(cs, 1.e-10f);
287}
288
289realnum AtomCSInterp(long int nelem,
290 long int ipHi,
291 long int ipLo,
292 realnum *factor1,
293 const char **where,
294 long int Collider )
295{
296 long ipArray;
297 realnum cs;
298
299 DEBUG_ENTRY( "AtomCSInterp()" );
300
301 ASSERT( nelem == ipHELIUM );
302
303 /* init values, better be positive when we exit */
304 cs = -1.f;
305
306 /* this may be used for splitting up the collision strength within 2 3P when
307 * the lower level is withint 2 3P, and for collisions between resolved and collapsed levels.
308 * It may be set somewhere in routine, so set to negative value here as flag saying not set */
309 *factor1 = -1.f;
310
311 /* for most of the helium iso sequence, the order of the J levels within 2 3P
312 * in increasing energy, is 0 1 2 - the exception is atomic helium itself,
313 * which is swapped, 2 1 0 */
314
315 /* this branch is for upper and lower levels within 2p3P */
316 if( ipLo >= ipHe2p3P0 && ipHi <= ipHe2p3P2 && Collider==ipELECTRON )
317 {
318 *factor1 = 1.f;
319 /* atomic helium, use Berrington private comm cs */
320
321 /* >>refer he1 cs Berrington, Keith, 2001, private communication - email follows
322 > Dear Gary,
323 > I could not find any literature on the He fine-structure
324 > problem (but I didn't look very hard, so there may be
325 > something somewhere). However, I did a quick R-matrix run to
326 > see what the magnitude of the collision strengths are... At
327 > 1000K, I get the effective collision strength for 2^3P J=0-1,
328 > 0-2, 1-2 as 0.8,0.7,2.7; for 10000K I get 1.2, 2.1, 6.0
329 */
330 /* indices are the same and correct, only thing special is that energies are in inverted order...was right first time. */
331 if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P1 )
332 {
333 cs = 1.2f;
334 }
335 else if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P2 )
336 {
337 cs = 2.1f;
338 }
339 else if( ipLo == ipHe2p3P1 && ipHi == ipHe2p3P2 )
340 {
341 cs = 6.0f;
342 }
343 else
344 {
345 cs = 1.0f;
347 }
348
349 *where = "Berr ";
350 /* statistical weights included */
351 }
352 /* >>chng 02 feb 25, Bray data should come first since it is the best we have. */
353 /* this branch is the Bray et al data, for n <= 5, where quantal calcs exist
354 * must exclude ipLo >= ipHe2p1P because they give no numbers for those */
355 else if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() <= 5 &&
356 ( ipHi < iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) &&
357 nelem==ipHELIUM && HeCS[ipHi][ipLo][0] >= 0.f && Collider== ipELECTRON )
358 {
359 ASSERT( *factor1 == -1.f );
360 ASSERT( ipLo < ipHi );
361 ASSERT( ipHe2p3P0 == 3 );
362
363 /* ipLo is within 2^3P */
364 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
365 {
366 /* *factor1 is ratio of statistical weights of level to term */
367
368 /* ipHe2p3P0, ipHe2p3P1, ipHe2p3P2 have indices 3,4,and 5, but j=0,1,and 2. */
369 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
370
371 /* ipHi must be above ipHe2p3P2 since 2p3Pj->2p3Pk taken care of above */
372 ASSERT( ipHi > ipHe2p3P2 );
373 }
374 /* ipHi is within 2^3P */
375 else if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
376 {
377 ASSERT( ipLo < ipHe2p3P0 );
378
379 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
380 }
381 /* neither are within 2^3P...no splitting necessary */
382 else
383 {
384 *factor1 = 1.f;
385 }
386
387 /* SOME OF THESE ARE NOT N-CHANGING! */
388 /* Must be careful about turning each one on or off. */
389
390 /* this is the case where we have quantal calculations */
391 /* >>refer He1 cs Bray, I., Burgess, A., Fursa, D.V., & Tully, J.A., 2000, A&AS, 146, 481-498 */
392 /* check whether we are outside temperature array bounds,
393 * and use extreme value if we are */
394 if( phycon.alogte <= CSTemp[0] )
395 {
396 cs = HeCS[ipHi][ipLo][0];
397 }
398 else if( phycon.alogte >= CSTemp.back() )
399 {
400 cs = HeCS[ipHi][ipLo][CSTemp.size()-1];
401 }
402 else
403 {
404 realnum flow;
405
406 /* find which array element within the cs vs temp array */
407 ipArray = (long)((phycon.alogte - CSTemp[0])/(CSTemp[1]-CSTemp[0]));
408 ASSERT( (unsigned)ipArray < CSTemp.size() );
409 ASSERT( ipArray >= 0 );
410 /* when taking the average, this is the fraction from the lower temperature value */
411 flow = (realnum)( (phycon.alogte - CSTemp[ipArray])/
412 (CSTemp[ipArray+1]-CSTemp[ipArray]));
413 ASSERT( (flow >= 0.f) && (flow <= 1.f) );
414
415 cs = HeCS[ipHi][ipLo][ipArray] * (1.f-flow) +
416 HeCS[ipHi][ipLo][ipArray+1] * flow;
417 }
418
419 *where = "Bray ";
420
421 /* options to kill collisional excitation and/or l-mixing */
422 if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() )
423 /* iso_ctrl.lgColl_l_mixing turned off with atom he-like l-mixing collisions off command */
424 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
425 else
426 {
427 /* iso_ctrl.lgColl_excite turned off with atom he-like collisional excitation off command */
428 cs *= (realnum)iso_ctrl.lgColl_excite[ipHE_LIKE];
429 }
430
431 ASSERT( cs >= 0.f );
432 /* statistical weights included */
433 }
434 /* this branch, n-same, l-changing collision, and not case of He with quantal data */
435 else if( (iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() ) &&
436 (iso_sp[ipHE_LIKE][nelem].st[ipHi].S() == iso_sp[ipHE_LIKE][nelem].st[ipLo].S() ) )
437 {
438 ASSERT( *factor1 == -1.f );
439 *factor1 = 1.f;
440
441 /* ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() >= 3 ); */
442 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max );
443
444 if( (iso_sp[ipHE_LIKE][nelem].st[ipLo].l() <=2) &&
445 abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1 )
446 {
447 /* Use the method given in
448 * >>refer He CS Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105
449 * statistical weights included */
450 cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
451 }
452 else if( iso_ctrl.lgCS_Vrinceanu[ipHE_LIKE] )
453 {
454 if( iso_sp[ipHE_LIKE][nelem].st[ipLo].l() >=3 &&
455 iso_sp[ipHE_LIKE][nelem].st[ipHi].l() >=3 )
456 {
457 /* Use the method given in
458 * >>refer He CS Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701
459 * statistical weights included */
461 nelem,
462 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
463 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
464 iso_sp[ipHE_LIKE][nelem].st[ipHi].l(),
465 iso_sp[ipHE_LIKE][nelem].st[ipLo].S(),
466 (double)phycon.te,
467 Collider );
468 }
469 else
470 {
471 cs = 0.f;
472 }
473 }
474 /* this branch, l changing by one */
475 else if( abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1)
476 {
477 /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
478 /* statistical weights included */
480 nelem,
481 iso_sp[ipHE_LIKE][nelem].st[ipLo].lifetime(),
482 nelem+1.-ipHE_LIKE,
483 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
484 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
485 iso_sp[ipHE_LIKE][nelem].st[ipHi].g(),
486 Collider);
487 }
488 else
489 {
490 /* l changes by more than 1, but same-n collision */
491 cs = 0.f;
492 }
493
494 /* ipLo is within 2^3P */
495 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
496 {
497 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
498 }
499
500 /* ipHi is within 2^3P */
501 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
502 {
503 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
504 }
505
506 *where = "lmix ";
507 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
508 }
509
510 /* this is an atomic n-changing collision with no quantal calculations */
511 else if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() != iso_sp[ipHE_LIKE][nelem].st[ipLo].n() )
512 {
513 ASSERT( *factor1 == -1.f );
514 /* this is an atomic n-changing collision with no quantal calculations */
515 /* gbar g-bar goes here */
516
517 /* >>chng 02 jul 25, add option for fits to quantal cs data */
518 if( iso_ctrl.lgCS_Vriens[ipHE_LIKE] )
519 {
520 /* >>refer He CS Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
521 * statistical weight IS included in the routine */
522 cs = (realnum)CS_VS80(
523 ipHE_LIKE,
524 nelem,
525 ipHi,
526 ipLo,
527 iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul(),
528 phycon.te,
529 Collider );
530
531 *factor1 = 1.f;
532 *where = "Vriens";
533 }
534 else if( iso_ctrl.lgCS_None[ipHE_LIKE] )
535 {
536 cs = 0.f;
537 *factor1 = 1.f;
538 *where = "no gb";
539 }
540 else if( iso_ctrl.nCS_new[ipHE_LIKE] )
541 {
542 *factor1 = 1.f;
543 /* Don't know if stat weights are included in this, but they're probably
544 * wrong anyway since they are based in part on the former (incorrect)
545 * implementation of Vriens and Smeets rates */
546
547 /* two different fits, allowed and forbidden */
548 if( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul() > 1. )
549 {
550 /* permitted lines - large A */
551 double x =
552 log10(MAX2(34.7,iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN()));
553
554 if( iso_ctrl.nCS_new[ipHE_LIKE] == 1 )
555 {
556 /* this is the broken power law fit, passing through both quantal
557 * calcs at high energy and asymptotically goes to VS at low energies */
558 if( x < 4.5 )
559 {
560 /* low energy fit for permitted transitions */
561 cs = (realnum)pow( 10. , -1.45*x + 6.75);
562 }
563 else
564 {
565 /* higher energy fit for permitted transitions */
566 cs = (realnum)pow( 10. , -3.33*x+15.15);
567 }
568 }
569 else if( iso_ctrl.nCS_new[ipHE_LIKE] == 2 )
570 {
571 /* single parallel fit for permitted transitions, runs parallel to VS */
572 cs = (realnum)pow( 10. , -2.3*x+10.3);
573 }
574 else
576 }
577 else
578 {
579 /* fit for forbidden transitions */
580 if( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN() < 25119.f )
581 {
582 cs = 0.631f;
583 }
584 else
585 {
586 cs = (realnum)pow(10.,
587 -3.*log10(iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN())+12.8);
588 }
589 }
590
591 *where = "newgb";
592 }
593 else
595
596 /* ipLo is within 2^3P */
597 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
598 {
599 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
600 }
601
602 /* ipHi is within 2^3P */
603 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
604 {
605 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
606 }
607
608 /* options to turn off collisions */
609 cs *= (realnum)iso_ctrl.lgColl_excite[ipHE_LIKE];
610
611 }
612 else
613 {
614 /* If spin changing collisions are prohibited in the l-mixing routine,
615 * they will fall here, and will have been assigned no collision strength.
616 * Assign zero for now. */
617 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].S() != iso_sp[ipHE_LIKE][nelem].st[ipLo].S() );
618 cs = 0.f;
619 *factor1 = 1.f;
620 }
621
622 ASSERT( cs >= 0.f );
623
624 return(cs);
625
626}
627
628/* IonCSInterp interpolate on collision strengths for element nelem */
629realnum IonCSInterp( long nelem , long ipHi , long ipLo, realnum *factor1, const char **where, long Collider )
630{
631 realnum cs;
632
633 DEBUG_ENTRY( "IonCSInterp()" );
634
635 ASSERT( nelem > ipHELIUM );
636 ASSERT( nelem < LIMELM );
637
638 /* init values, better be positive when we exit */
639 cs = -1.f;
640
641 /* this may be used for splitting up the collision strength for collisions between
642 * resolved and collapsed levels. It may be set somewhere in routine, so set to
643 * negative value here as flag saying not set */
644 *factor1 = -1.f;
645
646
647 /* >>chng 02 mar 04, the approximation here for transitions within 2p3P was not needed,
648 * because the Zhang data give these transitions. They are of the same order, but are
649 * specific to the three transitions */
650
651 /* this branch is ground to n=2 or from n=2 to n=2, for ions only */
652 /*>>refer Helike CS Zhang, Honglin, & Sampson, Douglas H. 1987, ApJS 63, 487 */
653 if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n()==2
654 && iso_sp[ipHE_LIKE][nelem].st[ipLo].n()<=2 && Collider==ipELECTRON)
655 {
656 *where = "Zhang";
657 *factor1 = 1.;
658
659 /* Collisions from gound */
660 if( ipLo == ipHe1s1S )
661 {
662 switch( ipHi )
663 {
664 case 1: /* to 2tripS */
665 cs = 0.25f/POW2(nelem+1.f);
666 break;
667 case 2: /* to 2singS */
668 cs = 0.4f/POW2(nelem+1.f);
669 break;
670 case 3: /* to 2tripP0 */
671 cs = 0.15f/POW2(nelem+1.f);
672 break;
673 case 4: /* to 2tripP1 */
674 cs = 0.45f/POW2(nelem+1.f);
675 break;
676 case 5: /* to 2tripP2 */
677 cs = 0.75f/POW2(nelem+1.f);
678 break;
679 case 6: /* to 2singP */
680 cs = 1.3f/POW2(nelem+1.f);
681 break;
682 default:
684 break;
685 }
686 cs *= (realnum)iso_ctrl.lgColl_excite[ipHE_LIKE];
687 }
688 /* collisions from 2tripS to n=2 */
689 else if( ipLo == ipHe2s3S )
690 {
691 switch( ipHi )
692 {
693 case 2: /* to 2singS */
694 cs = 2.75f/POW2(nelem+1.f);
695 break;
696 case 3: /* to 2tripP0 */
697 cs = 60.f/POW2(nelem+1.f);
698 break;
699 case 4: /* to 2tripP1 */
700 cs = 180.f/POW2(nelem+1.f);
701 break;
702 case 5: /* to 2tripP2 */
703 cs = 300.f/POW2(nelem+1.f);
704 break;
705 case 6: /* to 2singP */
706 cs = 5.8f/POW2(nelem+1.f);
707 break;
708 default:
710 break;
711 }
712 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
713 }
714 /* collisions from 2singS to n=2 */
715 else if( ipLo == ipHe2s1S )
716 {
717 switch( ipHi )
718 {
719 case 3: /* to 2tripP0 */
720 cs = 0.56f/POW2(nelem+1.f);
721 break;
722 case 4: /* to 2tripP1 */
723 cs = 1.74f/POW2(nelem+1.f);
724 break;
725 case 5: /* to 2tripP2 */
726 cs = 2.81f/POW2(nelem+1.f);
727 break;
728 case 6: /* to 2singP */
729 cs = 190.f/POW2(nelem+1.f);
730 break;
731 default:
733 break;
734 }
735 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
736 }
737 /* collisions from 2tripP0 to n=2 */
738 else if( ipLo == ipHe2p3P0 )
739 {
740 switch( ipHi )
741 {
742 case 4: /* to 2tripP1 */
743 cs = 8.1f/POW2(nelem+1.f);
744 break;
745 case 5: /* to 2tripP2 */
746 cs = 8.2f/POW2(nelem+1.f);
747 break;
748 case 6: /* to 2singP */
749 cs = 3.9f/POW2(nelem+1.f);
750 break;
751 default:
753 break;
754 }
755 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
756 }
757 /* collisions from 2tripP1 to n=2 */
758 else if( ipLo == ipHe2p3P1 )
759 {
760 switch( ipHi )
761 {
762 case 5: /* to 2tripP2 */
763 cs = 30.f/POW2(nelem+1.f);
764 break;
765 case 6: /* to 2singP */
766 cs = 11.7f/POW2(nelem+1.f);
767 break;
768 default:
770 break;
771 }
772 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
773 }
774 /* collisions from 2tripP2 to n=2 */
775 else if( ipLo == ipHe2p3P2 )
776 {
777 /* to 2singP */
778 cs = 19.5f/POW2(nelem+1.f) * (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
779 }
780 else
782
783 /* statistical weights included */
784 }
785
786 /* this branch, n-same, l-changing collisions */
787 else if( (iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() ) &&
788 (iso_sp[ipHE_LIKE][nelem].st[ipHi].S() == iso_sp[ipHE_LIKE][nelem].st[ipLo].S() ) )
789 {
790 ASSERT( *factor1 == -1.f );
791 *factor1 = 1.f;
792
793 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max );
794
795 if( (iso_sp[ipHE_LIKE][nelem].st[ipLo].l() <=2) &&
796 abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1 )
797 {
798 /* Use the method given in
799 * >>refer He CS Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105
800 * statistical weights included */
801 cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
802 }
803 else if( iso_ctrl.lgCS_Vrinceanu[ipHE_LIKE] )
804 {
805 if( iso_sp[ipHE_LIKE][nelem].st[ipLo].l() >=3 &&
806 iso_sp[ipHE_LIKE][nelem].st[ipHi].l() >=3 )
807 {
808 /* Use the method given in
809 * >>refer He CS Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701
810 * statistical weights included */
812 nelem,
813 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
814 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
815 iso_sp[ipHE_LIKE][nelem].st[ipHi].l(),
816 iso_sp[ipHE_LIKE][nelem].st[ipLo].S(),
817 (double)phycon.te,
818 Collider );
819 }
820 else
821 {
822 cs = 0.f;
823 }
824 }
825 /* this branch, l changing by one */
826 else if( abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1)
827 {
828 /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
829 /* statistical weights included */
831 nelem,
832 iso_sp[ipHE_LIKE][nelem].st[ipLo].lifetime(),
833 nelem+1.-ipHE_LIKE,
834 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
835 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
836 iso_sp[ipHE_LIKE][nelem].st[ipHi].g(),
837 Collider);
838 }
839 else
840 {
841 /* l changes by more than 1, but same-n collision */
842 cs = 0.f;
843 }
844
845 /* ipHi is within 2^3P, do not need to split on ipLo. */
846 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
847 {
848 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
849 }
850
851 *where = "lmix ";
852 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
853 }
854
855 /* this branch, n changing collisions for ions */
856 else if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() != iso_sp[ipHE_LIKE][nelem].st[ipLo].n() )
857 {
858 if( iso_ctrl.lgCS_Vriens[ipHE_LIKE] )
859 {
860 /* this is the default branch */
861 /* >>refer He CS Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
862 * statistical weight is NOT included in the routine */
863 cs = (realnum)CS_VS80(
864 ipHE_LIKE,
865 nelem,
866 ipHi,
867 ipLo,
868 iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul(),
869 phycon.te,
870 Collider );
871
872 *factor1 = 1.f;
873 *where = "Vriens";
874 }
875 else
876 {
877 /* \todo 2 this branch and above now do the same thing. Change something. */
878 /* this branch is for testing and reached with command ATOM HE-LIKE COLLISIONS VRIENS OFF */
879 fixit(); /* use Percival and Richards here */
880
881 cs = 0.f;
882 *where = "hydro";
883 }
884 }
885 else
886 {
887 /* what's left are deltaN=0, spin changing collisions.
888 * These have not been accounted for. */
889 /* Make sure what comes here is what we think it is. */
890 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() );
891 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].S() != iso_sp[ipHE_LIKE][nelem].st[ipLo].S() );
892 cs = 0.f;
893 *where = "spin ";
894 }
895
896 ASSERT( cs >= 0.f );
897
898 return(cs);
899}
900
901
902/*CS_l_mixing_S62 - find rate for l-mixing collisions by protons, for neutrals */
903/* The S62 stands for Seaton 1962 */
905 long ipISO,
906 long nelem /* the chemical element, 1 for He */,
907 long ipLo /* lower level, 0 for ground */,
908 long ipHi /* upper level, 0 for ground */,
909 double temp,
910 long Collider)
911{
912 /* >>refer He l-mixing Seaton, M.J., 1962, Proc. Phys. Soc. */
913 double coll_str;
914
915 DEBUG_ENTRY( "CS_l_mixing_S62()" );
916
917 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
918 {
919 return 0.;
920 }
921
922 my_Integrand_S62 func;
923
924 func.temp = temp;
925 func.stat_weight = iso_sp[ipISO][nelem].st[ipLo].g();
926 /* >>chng 05 sep 06, RP - update energies of excited states */
927 /* func.deltaE = EVRYD*(iso_sp[ipISO][nelem].st[ipLo].xIsoLevNIonRyd -
928 iso_sp[ipISO][nelem].st[ipHi].xIsoLevNIonRyd); */
929 func.deltaE = iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg()/EN1EV;
930 func.I_energy_eV = EVRYD*iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd;
931 func.Collider = Collider;
932 func.nelem = nelem;
933
935
936 func.osc_strength = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul()/
938
940 /* This returns a thermally averaged collision strength */
941 coll_str = S62.sum( 0., 1., func );
942 coll_str += S62.sum( 1., 10., func );
943 ASSERT( coll_str > 0. );
944
945 return coll_str;
946}
947
948/* The integrand for calculating the thermal average of collision strengths */
949STATIC double S62_Therm_ave_coll_str( double proj_energy_overKT, long nelem, long Collider, double deltaE, double osc_strength, double temp,
950 double stat_weight, double I_energy_eV )
951{
952
953 double integrand, cross_section, /*Rnot,*/ coll_str, zOverB2;
954 double /* betanot, */ betaone, zeta_of_betaone, /* cs1, */ cs2;
955 double Dubya, proj_energy;
956 double reduced_mass;
957
958 DEBUG_ENTRY( "S62_Therm_ave_coll_str()" );
959
960 /* projectile energy in eV */
961 proj_energy = proj_energy_overKT * phycon.te / EVDEGK;
962
963 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
964 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
965
966 /* Rnot = 1.1229*H_BAR/sqrt(ELECTRON_MASS*deltaE*EN1EV)/Bohr_radius; in units of Bohr_radius */
967
968 proj_energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
969 /* The deltaE here is to make sure that the collider has no less
970 * than the energy difference between the initial and final levels. */
971 proj_energy += deltaE;
972 Dubya = 0.5*(2.*proj_energy-deltaE);
973 ASSERT( Dubya > 0. );
974
975 /* betanot = sqrt(proj_energy/I_energy_eV)*(deltaE/2./Dubya)*Rnot; */
976
977 ASSERT( I_energy_eV > 0. );
978 ASSERT( osc_strength > 0. );
979
980 /* from equation 33 */
981 zOverB2 = 0.5*POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength;
982
983 ASSERT( zOverB2 > 0. );
984
985 if( zOverB2 > 100. )
986 {
987 betaone = sqrt( 1./zOverB2 );
988 }
989 else if( zOverB2 < 0.54 )
990 {
991 /* Low betaone approximation */
992 betaone = (1./3.)*(log(PI)-log(zOverB2)+1.28);
993 if( betaone > 2.38 )
994 {
995 /* average with this over approximation */
996 betaone += 0.5*(log(PI)-log(zOverB2));
997 betaone *= 0.5;
998 }
999 }
1000 else
1001 {
1002 long ip_zOverB2 = 0;
1003 double zetaOVERbeta2[101] = {
1004 1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
1005 7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
1006 4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
1007 3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
1008 2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
1009 1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
1010 1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
1011 7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
1012 4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
1013 3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
1014 2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
1015 1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
1016 7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};
1017
1018 ASSERT( zOverB2 >= zetaOVERbeta2[100] );
1019
1020 /* find beta in the table */
1021 for( long i=0; i< 100; ++i )
1022 {
1023 if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) )
1024 {
1025 /* found the temperature - use it */
1026 ip_zOverB2 = i;
1027 break;
1028 }
1029 }
1030
1031 ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );
1032
1033 betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) /
1034 (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) *
1035 (pow(10., ((double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((double)ip_zOverB2)/100. - 1.)) +
1036 pow(10., (double)ip_zOverB2/100. - 1.);
1037
1038 }
1039
1040 zeta_of_betaone = zOverB2 * POW2(betaone);
1041
1042 /* cs1 = betanot * bessel_k0(betanot) * bessel_k1(betanot); */
1043 cs2 = 0.5*zeta_of_betaone + betaone * bessel_k0(betaone) * bessel_k1(betaone);
1044
1045 /* cross_section = MIN2(cs1, cs2); */
1046 cross_section = cs2;
1047
1048 /* cross section in units of PI * a_o^2 */
1049 cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy);
1050
1051 /* convert to collision strength */
1052 coll_str = ConvCrossSect2CollStr( cross_section * PI*BOHR_RADIUS_CM*BOHR_RADIUS_CM, stat_weight, proj_energy/EVRYD, reduced_mass );
1053
1054 integrand = exp( -1.*(proj_energy-deltaE)*EVDEGK/temp ) * coll_str;
1055 return integrand;
1056}
1057
1058/*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
1060 long nelem, /* the chemical element, 1 for He */
1061 double tau, /* the radiative lifetime of the initial level. */
1062 double target_charge,
1063 long n,
1064 long l,
1065 double gHi,
1066 long Collider)
1067{
1068 /* >>refer H-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
1069 /* >>refer He-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
1070 double cs, Dnl,
1071 TwoLogDebye, TwoLogRc1,
1072 factor1, factor2,
1073 bestfactor, factorpart,
1074 reduced_mass, reduced_mass_2_emass,
1075 rate;
1076 const double ChargIncoming = ColliderCharge[Collider];
1077
1078 DEBUG_ENTRY( "CS_l_mixing_PS64()" );
1079
1080 /* In this routine, two different cutoff radii are calculated, and as per PS64,
1081 * the least of these is selected. We take the least positive result. */
1082
1083 /* This reduced mass is in grams. */
1084 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
1085 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
1086 /* this mass always appears relative to the electron mass, so define it that way */
1087 reduced_mass_2_emass = reduced_mass / ELECTRON_MASS;
1088
1089 /* equation 46 of PS64 */
1090 /* min on density added to prevent this from becoming large and negative
1091 * at very high n_e - Pengelly & Seaton did not apply this above 1e12 cm^-3 */
1092 /* This is 2 times the log of the Debye radius. */
1093 TwoLogDebye = 1.68 + log10( phycon.te / MIN2(1e11 , dense.eden ) );
1094 /* Brocklehurst (1971, equation 3.40) has 1.181 instead of 1.68. This appears to be due to
1095 * Pengelly and Seaton neglecting one of the two factors of PI in their Equation 33 */
1096 //TwoLogDebye = 1.181 + log10( phycon.te / MIN2(1e10 , dense.eden ) );
1097
1098 /* This is 2 times the log of cutoff = 0.72v(tau), where tau is the lifetime of the level nl.
1099 * This is PS64 equation 45 (same as Brocklehurst 1971 equation 3.41) */
1100 TwoLogRc1 = 10.95 + log10( phycon.te * tau * tau / reduced_mass_2_emass );
1101
1102#if 0
1103 /* This is 2 times the log of cutoff = 1.12 * hbar v / deltaE, where v is reduced velocity = sqrt( 2kT/mu ). */
1104 /* This is PS64 equation 29 */
1105 if( deltaE > 0. )
1106 TwoLogRc2 = 2. * log10( 1.12 * H_BAR * v / deltaE );
1107 else
1108 TwoLogRc2 = BIGDOUBLE;
1109#endif
1110
1111 /* this is equation 44 of PS64 */
1112 Dnl = POW2( ChargIncoming / target_charge) * 6. * POW2( (double)n) *
1113 ( POW2((double)n) - POW2((double)l) - l - 1);
1114
1115 ASSERT( Dnl > 0. );
1116 ASSERT( phycon.te / Dnl / reduced_mass_2_emass > 0. );
1117
1118 factorpart = (11.54 + log10( phycon.te / Dnl / reduced_mass_2_emass ) );
1119
1120 if( (factor1 = factorpart + TwoLogDebye) <= 0.)
1121 factor1 = BIGDOUBLE;
1122 if( (factor2 = factorpart + TwoLogRc1) <= 0.)
1123 factor2 = BIGDOUBLE;
1124
1125 /* Now we find the least positive result. */
1126 bestfactor = MIN2(factor1,factor2);
1127
1128 ASSERT( bestfactor > 0. );
1129
1130 /* If both factors are bigger than 100, toss out the result, and return SMALLFLOAT. */
1131 if( bestfactor > 100. )
1132 {
1133 return SMALLFLOAT;
1134 }
1135
1136 /* This is the rate coefficient. Units: cm^3 s-1 */
1137 rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dnl / phycon.sqrte * bestfactor;
1138
1139 /***** NB NB NB NB
1140 * Brocklehurst (1971) has a factor of electron density in the rate coefficient (equation 3.38).
1141 * This is NOT a proper rate, as can be seen in his equations 2.2 and 2.4. This differs
1142 * from the formulism given by PS64. */
1143 //rate *= dense.eden;
1144
1145 /* this is the TOTAL rate from nl to nl+/-1. So assume we can
1146 * divide by two to get the rate nl to either nl+1 or nl-1. */
1147 if( l > 0 )
1148 rate /=2;
1149
1150 /* convert rate to collision strength */
1151 /* NB - the term in parentheses corrects for the fact that COLL_CONST is only appropriate
1152 * for electron colliders and is off by reduced_mass_2_emass^-1.5 */
1153 cs = rate / ( COLL_CONST * pow(reduced_mass_2_emass, -1.5) ) *
1154 phycon.sqrte * gHi;
1155
1156 ASSERT( cs > 0. );
1157
1158 return cs;
1159}
1160
1161/*CS_l_mixing - find collision strength for l-mixing collisions for neutrals */
1162/* The VF stands for Vrinceanu & Flannery 2001 */
1163/* >>refer He l-mixing Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701 */
1164/* >>refer He l-mixing Hezel, T. P., Burkhardt, C. E., Ciocca, M., He, L-W., */
1165/* >>refercon Leventhal, J. J. 1992, Am. J. Phys. 60, 329 */
1166double CS_l_mixing_VF01(long int ipISO,
1167 long int nelem,
1168 long int n,
1169 long int l,
1170 long int lp,
1171 long int s,
1172 double temp,
1173 long int Collider )
1174{
1175
1176 double coll_str;
1177
1178 DEBUG_ENTRY( "CS_l_mixing_VF01()" );
1179
1181 func.ipISO = ipISO;
1182 func.nelem = nelem;
1183 func.n = n;
1184 func.l = l;
1185 func.lp = lp;
1186 func.s = s;
1187 func.Collider = Collider;
1188 func.temp = temp;
1189 func.ColliderCharge = ColliderCharge[Collider];
1190 func.lgParamIsRedVel = false;
1191 ASSERT( func.ColliderCharge > 0. );
1192
1194
1195 /* no need to do this for h-like */
1196 if( ipISO > ipH_LIKE )
1197 {
1198 ASSERT( l != 0 );
1199 ASSERT( lp != 0 );
1200 }
1201
1202 if( !iso_ctrl.lgCS_therm_ave[ipISO] )
1203 {
1204 /* Must do some thermal averaging for densities greater
1205 * than about 10000 and less than about 1e10,
1206 * because kT gives significantly different results.
1207 * Still, do sparser integration than is done below */
1208 if( (dense.eden > 10000.) && (dense.eden < 1E10 ) )
1209 {
1210 coll_str = VF01_E.sum( 0.0, 6.0, func );
1211 }
1212 else
1213 {
1214 /* Do NOT average over Maxwellian */
1215 coll_str = collision_strength_VF01( ipISO, nelem, n, l, lp, s, Collider,
1216 ColliderCharge[Collider], temp, temp/TE1RYD, false );
1217 }
1218 }
1219 else
1220 {
1221 /* DO average over Maxwellian */
1222 coll_str = VF01_E.sum( 0.0, 1.0, func );
1223 coll_str += VF01_E.sum( 1.0, 10.0, func );
1224 }
1225
1226 return coll_str;
1227}
1228
1229STATIC double collision_strength_VF01( long ipISO, long nelem, long n, long l, long lp, long s, long Collider, double ColliderCharge, double temp, double velOrEner, bool lgParamIsRedVel )
1230{
1231 double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd;
1232 double ConstantFactors, reduced_mass, CSIntegral, stat_weight;
1233 double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd;
1234 double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 /*, alpha2*/;
1235 long ipLo, ipHi;
1236
1237 DEBUG_ENTRY( "collision_strength_VF01()" );
1238
1239 ASSERT( n > 0 );
1240
1241 /* >>chng 06 may 30, move these up from below. Also ipHi needs lp not l. */
1242 ipLo = iso_sp[ipISO][nelem].QuantumNumbers2Index[n][l][s];
1243 ipHi = iso_sp[ipISO][nelem].QuantumNumbers2Index[n][lp][s];
1244 stat_weight = iso_sp[ipISO][nelem].st[ipLo].g();
1245
1246 /* no need to do this for h-like */
1247 if( ipISO > ipH_LIKE )
1248 {
1249 /* these shut up the lint, already done above */
1250 ASSERT( l > 0 );
1251 ASSERT( lp > 0 );
1252 }
1253
1254 /* This reduced mass is in grams. */
1255 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
1256 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
1257 ASSERT( reduced_mass > 0. );
1258
1259 /* this is root mean squared velocity */
1260 /* use this as projectile velocity for thermally-averaged cross-section */
1261 aveRadius = (BOHR_RADIUS_CM/((double)nelem+1.-(double)ipISO))*POW2((double)n);
1262 ASSERT( aveRadius < 1.e-4 );
1263 /* >>chng 05 jul 14, as per exchange with Ryan Porter & Peter van Hoof, avoid
1264 * roundoff error and give ability to go beyond zinc */
1265 /*ASSERT( aveRadius >= BOHR_RADIUS_CM );*/
1266 ASSERT( aveRadius > 3.9/LIMELM * BOHR_RADIUS_CM );
1267
1268 /* vn = n * H_BAR/ m / r = Z * e^2 / n / H_BAR
1269 * where Z is the effective charge. */
1270 RMSv = ((double)nelem+1.-(double)ipISO)*POW2(ELEM_CHARGE_ESU)/(double)n/H_BAR;
1271 ASSERT( RMSv > 0. );
1272
1273 ASSERT( ColliderMass[Collider] > 0. );
1274
1275 if( lgParamIsRedVel )
1276 {
1277 /* velOrEner is a reduced velocity */
1278 reduced_vel = velOrEner;
1279 /* The proton mass is needed here because the ColliderMass array is
1280 * expressed in units of the proton mass, but here we need absolute mass. */
1281 E_Proj_Ryd = 0.5 * POW2( reduced_vel * RMSv ) * ColliderMass[Collider] *
1283 }
1284 else
1285 {
1286 /* velOrEner is a projectile energy in Rydbergs */
1287 E_Proj_Ryd = velOrEner;
1288 reduced_vel = sqrt( 2.*E_Proj_Ryd*EN1RYD/ColliderMass[Collider]/PROTON_MASS )/RMSv;
1289 }
1290
1291 /* put limits on the reduced velocity. These limits should be more than fair. */
1292 ASSERT( reduced_vel > 1.e-10 );
1293 ASSERT( reduced_vel < 1.e10 );
1294
1295 /* Factors outside integral */
1296 ConstantFactors = 4.5*PI*POW2(ColliderCharge*aveRadius/reduced_vel);
1297
1298 /* Reduced here means in units of aveRadius: */
1299 reduced_b_min = 1.5 * ColliderCharge / reduced_vel;
1300 ASSERT( reduced_b_min > 1.e-10 );
1301 ASSERT( reduced_b_min < 1.e10 );
1302
1303 if( ipISO == ipH_LIKE )
1304 {
1305 /* Debye radius: appears to be too large, results in 1/v^2 variation. */
1306 reduced_b_max = sqrt( BOLTZMANN*temp/ColliderCharge/dense.eden )
1307 / (PI2*ELEM_CHARGE_ESU)/aveRadius;
1308 }
1309 else if( ipISO == ipHE_LIKE )
1310 {
1311 quantum_defect1 = (double)n- (double)nelem/sqrt(iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd);
1312 quantum_defect2 = (double)n- (double)nelem/sqrt(iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
1313
1314 /* The magnitude of each quantum defect must be between zero and one. */
1315 ASSERT( fabs(quantum_defect1) < 1.0 );
1316 ASSERT( fabs(quantum_defect1) > 0.0 );
1317 ASSERT( fabs(quantum_defect2) < 1.0 );
1318 ASSERT( fabs(quantum_defect2) > 0.0 );
1319
1320 /* The quantum defect precession frequencies */
1321 omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l );
1322 /* >>chng 06 may 30, this needs lp not l. */
1323 omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*POW2((double)lp/(double)n)) / POW3( (double)n ) / (double)lp );
1324 /* Take the average for the two levels, for reciprocity. */
1325 omega_qd = 0.5*( omega_qd1 + omega_qd2 );
1326
1327 ASSERT( omega_qd > 0. );
1328
1329 reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius;
1330 }
1331 else
1332 /* rethink this before using on other iso sequences. */
1333 TotalInsanity();
1334
1335 reduced_b_max = MAX2( reduced_b_max, reduced_b_min );
1336 ASSERT( reduced_b_max > 0. );
1337
1338 // set up the integrator.
1340 func.n = n;
1341 func.l = l;
1342 func.lp = lp;
1343 func.red_vel = reduced_vel;
1344 func.an = aveRadius;
1346 func.bmax = reduced_b_max*aveRadius;
1348
1349 alphamin = 1.5*ColliderCharge/(reduced_vel * reduced_b_max);
1350 alphamax = 1.5*ColliderCharge/(reduced_vel * reduced_b_min);
1351
1352 ASSERT( alphamin > 0. );
1353 ASSERT( alphamax > 0. );
1354
1355 alphamin = MAX2( alphamin, 1.e-30 );
1356 alphamax = MAX2( alphamax, 1.e-20 );
1357
1358 CSIntegral = 0.;
1359
1360 if( alphamax > alphamin )
1361 {
1362
1363 step = (alphamax - alphamin)/5.;
1364 alpha1 = alphamin;
1365 CSIntegral += VF01_alpha.sum( alpha1, alpha1+step, func );
1366 CSIntegral += VF01_alpha.sum( alpha1+step, alpha1+4.*step, func );
1367 }
1368
1369 /* Calculate cross section */
1370 cross_section = ConstantFactors * CSIntegral;
1371
1372 /* convert to collision strength, cross section is already in cm^2 */
1373 coll_str = ConvCrossSect2CollStr( cross_section, stat_weight, E_Proj_Ryd, reduced_mass );
1374
1375 coll_str = MAX2( SMALLFLOAT, coll_str);
1376
1377 return coll_str;
1378}
1379
1380STATIC double L_mix_integrand_VF01( long n, long l, long lp, double bmax, double red_vel, double an, double ColliderCharge, double alpha )
1381{
1382 double integrand, deltaPhi, b;
1383
1384 DEBUG_ENTRY( "L_mix_integrand_VF01()" );
1385
1386 ASSERT( alpha >= 1.e-30 );
1387 ASSERT( bmax > 0. );
1388 ASSERT( red_vel > 0. );
1389
1390 /* >>refer He l-mixing Kazansky, A. K. & Ostrovsky, V. N. 1996, JPhysB: At. Mol. Opt. Phys. 29, 3651 */
1391 b = 1.5*ColliderCharge*an/(red_vel * alpha);
1392 /* deltaPhi is the effective angle swept out by the projector as viewed by the target. */
1393 if( b < bmax )
1394 {
1395 deltaPhi = -1.*PI + 2.*asin(b/bmax);
1396 }
1397 else
1398 {
1399 deltaPhi = 0.;
1400 }
1401 integrand = 1./alpha/alpha/alpha;
1402 integrand *= StarkCollTransProb_VF01( n, l, lp, alpha, deltaPhi );
1403
1404 return integrand;
1405}
1406
1407STATIC double StarkCollTransProb_VF01( long n, long l, long lp, double alpha, double deltaPhi)
1408{
1409 double probability;
1410 double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B;
1411
1412 DEBUG_ENTRY( "StarkCollTransProb_VF01()" );
1413
1414 ASSERT( alpha > 0. );
1415
1416 /* These are defined on page 11 of VF01 */
1417 cosU1 = 2.*POW2((double)l/(double)n) - 1.;
1418 cosU2 = 2.*POW2((double)lp/(double)n) - 1.;
1419
1420 sinU1 = sqrt( 1. - cosU1*cosU1 );
1421 sinU2 = sqrt( 1. - cosU2*cosU2 );
1422
1423
1424 cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha);
1425 sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 );
1426 cosChi = 2. * POW2( cosChiOver2 ) - 1.;
1427
1428 if( l == 0 )
1429 {
1430 if( -1.*cosU2 - cosChi < 0. )
1431 {
1432 probability = 0.;
1433 }
1434 else
1435 {
1436 /* Here is the initial state S case */
1437 ASSERT( sinChiOver2 > 0. );
1438 ASSERT( sinChiOver2*sinChiOver2 > POW2((double)lp/(double)n) );
1439 /* This is equation 35 of VF01. There is a factor of hbar missing in the denominator of the
1440 * paper, but it's okay if you use atomic units (hbar = 1). */
1441 probability = (double)lp/(POW2((double)n)*sinChiOver2*sqrt( POW2(sinChiOver2) - POW2((double)lp/(double)n) ) );
1442 }
1443 }
1444 else
1445 {
1446 double OneMinusCosChi = 1. - cosChi;
1447
1448 if( OneMinusCosChi == 0. )
1449 {
1450 double hold = sin( deltaPhi / 2. );
1451 /* From approximation at bottom of page 10 of VF01. */
1452 OneMinusCosChi = 8. * alpha * alpha * POW2( hold );
1453 }
1454
1455 if( OneMinusCosChi == 0. )
1456 {
1457 probability = 0.;
1458 }
1459 else
1460 {
1461 /* Here is the general case */
1462 A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi;
1463 B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi;
1464
1465 ASSERT( B > A );
1466
1467 /* These are the three cases of equation 34. */
1468 if( B <= 0. )
1469 {
1470 probability = 0.;
1471 }
1472 else
1473 {
1474 ASSERT( POW2( sinChiOver2 ) > 0. );
1475
1476 probability = 2.*lp/(PI* /*H_BAR* */ n*n*POW2( sinChiOver2 ));
1477
1478 if( A < 0. )
1479 {
1480 probability *= ellpk( -A/(B-A) );
1481 probability /= sqrt( B - A );
1482 }
1483 else
1484 {
1485 probability *= ellpk( A/B );
1486 probability /= sqrt( B );
1487 }
1488 }
1489 }
1490
1491 }
1492
1493 return probability;
1494
1495}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define POW3
Definition cddefines.h:936
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const char * strchr_s(const char *s, int c)
Definition cddefines.h:1439
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
double sum(double min, double max, Integrand func)
Definition cddefines.h:1531
double operator()(double proj_energy_overKT)
Definition helike_cs.cpp:48
double operator()(double EOverKT)
Definition helike_cs.cpp:63
double operator()(double alpha)
Definition helike_cs.cpp:77
@ ipELECTRON
Definition collision.h:9
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const double BIGDOUBLE
Definition cpu.h:194
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
#define COLLISMAGIC
Definition helike.h:24
STATIC double L_mix_integrand_VF01(long n, long l, long lp, double bmax, double red_vel, double an, double ColliderCharge, double alpha)
STATIC double S62_Therm_ave_coll_str(double proj_energy_overKT, long nelem, long Collider, double deltaE, double osc_strength, double temp, double stat_weight, double I_energy_eV)
STATIC double StarkCollTransProb_VF01(long int n, long int l, long int lp, double alpha, double deltaPhi)
realnum HeCSInterp(long int nelem, long int ipHi, long int ipLo, long int Collider)
vector< double > CSTemp
Definition helike_cs.cpp:26
realnum IonCSInterp(long nelem, long ipHi, long ipLo, realnum *factor1, const char **where, long Collider)
double CS_l_mixing_VF01(long int ipISO, long int nelem, long int n, long int l, long int lp, long int s, double temp, long int Collider)
static const double ColliderCharge[4]
Definition helike_cs.cpp:88
double CS_l_mixing_S62(long ipISO, long nelem, long ipLo, long ipHi, double temp, long Collider)
void HeCollidSetup(void)
Definition helike_cs.cpp:90
#define chLine_LENGTH
multi_arr< realnum, 3 > HeCS
Definition helike_cs.cpp:28
double CS_l_mixing_PS64(long nelem, double tau, double target_charge, long n, long l, double gHi, long Collider)
realnum AtomCSInterp(long int nelem, long int ipHi, long int ipLo, realnum *factor1, const char **where, long int Collider)
static const double ColliderMass[4]
Definition helike_cs.cpp:87
STATIC double collision_strength_VF01(long ipISO, long nelem, long n, long l, long lp, long s, long Collider, double ColliderCharge, double temp, double velOrEner, bool lgParamIsRedVel)
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
double CS_VS80(long int ipISO, long int nelem, long int ipHi, long int ipLo, double Aul, double temp, long int Collider)
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
const int ipHe2s3S
Definition iso.h:44
const int ipHe2s1S
Definition iso.h:45
const int ipHE_LIKE
Definition iso.h:63
const int ipHe2p3P1
Definition iso.h:47
const int ipHe2p3P0
Definition iso.h:46
const int ipH_LIKE
Definition iso.h:62
const int ipHe2p3P2
Definition iso.h:48
double ConvCrossSect2CollStr(double CrsSectCM2, double gLo, double E_ProjectileRyd, double reduced_mass_grams)
#define S(I_, J_)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double H_BAR
Definition physconst.h:144
UNUSED const double TRANS_PROB_CONST
Definition physconst.h:237
UNUSED const double PI
Definition physconst.h:29
UNUSED const double BOHR_RADIUS_CM
Definition physconst.h:222
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double PI2
Definition physconst.h:32
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double EN1EV
Definition physconst.h:192
UNUSED const double PROTON_MASS
Definition physconst.h:94
UNUSED const double EVDEGK
Definition physconst.h:186
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
UNUSED const double COLL_CONST
Definition physconst.h:229
UNUSED const double ELEM_CHARGE_ESU
Definition physconst.h:147
UNUSED const double WAVNRYD
Definition physconst.h:173
UNUSED const double TE1RYD
Definition physconst.h:183
static double ** col_str
Definition species2.cpp:29
static double * g
Definition species2.cpp:28
double bessel_k0(double x)
double ellpk(double x)
double bessel_k1(double x)
t_trace trace
Definition trace.cpp:5