45 long int i, nelem, ipHi, ipLo;
50 photons_3889_plus_7065 = 0.;
55 fprintf(
ioQQQ,
" prt_lines_helium called\n" );
63 " start He-like iso sequence");
66 " total collisional cooling due to all HeI lines ");
69 " total collisional heating due to all HeI lines ");
79 for( nelem=ipISO; nelem <
LIMELM; nelem++ )
81 if(
dense.lgElmtOn[nelem] )
94 for(
long ipDens = 0; ipDens <
NUMDENS; ++ipDens )
97 intens = pow( 10., intens ) *
dense.xIonDense[nelem][nelem+1-ipISO]*
dense.eden;
108 chIonLbl(chLabel, nelem+1, nelem+1-ipISO);
110 for( vector<two_photon>::iterator tnu = sp->
TwoNu.begin(); tnu != sp->
TwoNu.end(); ++tnu )
114 linadd( tnu->AulTotal * tnu->E2nu *
EN1RYD * (*tnu->Pop),
116 " two photon continuum ");
118 linadd( tnu->induc_dn * tnu->E2nu *
EN1RYD * (*tnu->Pop),
120 " induced two photon emission ");
141 " total emission in He-like intercombination lines from 2p3P to ground ");
147 if(
prt.lgPrnIsoCollapsed )
156 vector<long> EnterTheseLast;
157 for( ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
166 if(
iso_ctrl.lgFSM[ipISO] && ( abs(sp->
st[ipHi].l() -
167 sp->
st[ipLo].l())==1 )
168 && (sp->
st[ipLo].l()>1)
169 && (sp->
st[ipHi].l()>1)
170 && ( sp->
st[ipHi].n() ==
174 if( (sp->
st[ipHi].S()==1)
175 && (sp->
st[ipLo].S()==1) )
179 else if( (sp->
st[ipHi].S()==3)
180 && (sp->
st[ipLo].S()==3) )
187 (sp->
trans(ipHi,ipLo+1).Emis().Pesc() +
190 sp->
st[ipHi+1].Pop()*
191 (sp->
trans(ipHi+1,ipLo+1).Emis().Pesc() +
204 (sp->
trans(ipHi,ipLo).Emis().Pesc() +
207 sp->
st[ipHi+1].Pop()*
208 (sp->
trans(ipHi+1,ipLo).Emis().Pesc() +
233 (sp->
trans(i,ipLo).Emis().Pesc() +
243 fprintf(
ioQQQ,
"DEBUG 2P - 2S multiplet wl %s ",
246 fprintf(
ioQQQ,
"\n" );
250 linadd(sum,av_wl,
"TOTL",
'i',
251 "total emission in He-like lines, use average of three line wavelengths " );
254 linadd(sum,av_wl,chLabel,
'i',
255 "total emission in He-like lines, use average of three line wavelengths " );
264 (sp->
trans(ipHi,ipLo).Emis().Pesc() +
275 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
277 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
288 (sp->
trans(ipHi,ipLo).Emis().Pesc() +
300 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
302 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
305 if( abs(
L_(ipHi) -
L_(ipLo) ) != 1 )
307 EnterTheseLast.push_back( ipHi );
315 "total intensity of He-like line");
319 enum {DEBUG_LOC=
false};
322 if( nelem==1 && ipLo==0 && ipHi==1 )
324 fprintf(
ioQQQ,
"he1 626 %.2e %.2e \n",
335 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
337 "predicted line, all processes included");
341 for( ipHi=
ipHe2p3P2+1; ipHi < nLoop; ipHi++ )
343 double sumcool , sumheat ,
344 save , savecool , saveheat;
358 (sp->
trans(ipHi,ipLo).Emis().Pesc() +
370 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
372 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
396 "predicted line, all processes included");
402 for( ipLo=
ipHe2p3P2+1; ipLo < nLoop-1; ipLo++ )
404 vector<long> EnterTheseLast;
405 for( ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
415 (sp->
trans(ipHi,ipLo).Emis().Pesc() +
427 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
429 sp->
ex[ipHi][ipLo].ErrorFactor[
IPRAD];
432 if( abs(
L_(ipHi) -
L_(ipLo) ) != 1 )
434 EnterTheseLast.push_back( ipHi );
440 "predicted line, all processes included");
444 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
446 "predicted line, all processes included");
464 photons_3889_plus_7065 =
495 photons_3889_plus_7065 +=
503 linadd( photons_3889_plus_7065, 3889.,
"Pho+",
'i',
504 "photon sum given in Porter et al. 2007 (astro-ph/0611579)");
513 fprintf(
ioQQQ,
" lines_helium returns\n" );
525# define chLine_LENGTH 1000
531 fprintf(
ioQQQ,
" lines_helium opening he1_case_b.dat\n");
533 ioDATA =
open_data(
"he1_case_b.dat",
"r" );
538 fprintf(
ioQQQ,
" lines_helium could not read first line of he1_case_b.dat.\n");
543 i1 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
544 i2 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
551 " lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
553 " lines_helium: I expected to find the number %i and got %li instead.\n" ,
555 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
561 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
564 if( chLine[0] !=
'#')
581 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
584 if( chLine[0] !=
'#')
601 fprintf(
ioQQQ,
" lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
627 for(
long j = 0; j < i2; ++j )
633 strcpy(
CaBLines[nelem][j].label ,
" " );
635 for(
long k = 0; k <
NUMDENS; ++k )
638 for(
long l = 0; l <
NUMTEMPS; ++l )
651 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
656 if( (chLine[0] ==
' ') || (chLine[0]==
'\n') )
658 if( chLine[0] !=
'#')
667 long ipLo = (long)
FFmtRead(chLine,&j,
sizeof(chLine),&lgEOL);
668 long ipHi = (long)
FFmtRead(chLine,&j,
sizeof(chLine),&lgEOL);
674 strncpy(
CaBLines[nelem][line].label,
"Ca B" , 4 );
681 atmdat.CaseBWlHeI.push_back( wl );
683 for(
long ipDens = 0; ipDens <
NUMDENS; ++ipDens )
687 fprintf(
ioQQQ,
" lines_helium could not scan case B lines, current indices: %li %li\n",
695 long den = (long)
FFmtRead(chTemp,&j,
sizeof(chTemp),&lgEOL);
696 ASSERT( den == ipDens + 1 );
698 for(
long ipTe = 0; ipTe <
NUMTEMPS; ++ipTe )
701 if( (chTemp =
strchr_s( chTemp,
'\t' )) == NULL )
703 fprintf(
ioQQQ,
" lines_helium could not scan case B lines, current indices: %li %li\n",
709 sscanf( chTemp,
"%le" , &b );
735 if( Te <= TempArray[0] )
737 return ValueArray[0] + log10( TempArray[0] ) - log10( Te );
739 else if( Te >= TempArray[NumElements-1] )
741 return ValueArray[NumElements-1];
747 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
751 i0 =
max(
min(ipTe-ORDER/2,NumElements-ORDER-1),0);
752 rate =
lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, Te );