32 STATIC double TempInterp2(
double* TempArray ,
double* ValueArray,
long NumElements,
double Te );
48 trList.push_back( tr );
55 double &sum_inten,
double &sum_obs_inten,
double &sum_cool,
double &sum_heat )
58 sum_inten = sum_obs_inten = sum_cool = sum_heat = 0.;
60 if( trList.size() == 0)
64 for( vector<TransitionProxy>::iterator tr = trList.begin(); tr != trList.end(); tr++ )
66 sum_obs_inten += (*tr).Emis().xObsIntensity();
67 sum_inten += (*tr).Emis().xIntensity();
68 sum_cool += (*tr).Coll().cool();
69 sum_heat += (*tr).Coll().heat();
70 wl.push_back( (*tr).WLAng() );
76 std::sort( wl.begin(), wl.end(), std::less<realnum>() );
86 const double sumxInt,
const double sumxObsInt,
87 const double sumcool,
const double sumheat,
88 const char *chComment )
110 for( vector<TransitionProxy>::iterator tr = multiHe.begin(); tr != multiHe.end(); tr++ )
112 if( (*tr).ipHi() == ipHi && (*tr).ipLo() == ipLo )
148 " start He-like iso sequence");
158 for(
long nelem=ipISO; nelem <
LIMELM; nelem++ )
160 vector<TransitionProxy> multiplet;
161 vector<TransitionProxy> HeTrList;
162 double sumxInt = 0., sumxObsInt = 0., sumcool = 0., sumheat = 0.;
177 chLabel =
chIonLbl(nelem+1, nelem+1-ipISO);
179 for( vector<two_photon>::iterator tnu = sp->
TwoNu.begin(); tnu != sp->
TwoNu.end(); ++tnu )
181 fixit(
"This was multiplied by Pesc when treated as a line, now what? Only used for printout?");
182 fixit(
"below should be 'i' instead of 'r' ?");
184 string tpc_comment =
"",
189 tpc_comment =
" two photon continuum, " + comment_trans;
190 tpe_comment =
" induced two photon emission, " + comment_trans;
193 linadd( tnu->AulTotal * tnu->E2nu * EN1RYD * (*tnu->Pop),
194 0, chLabel.c_str(),
'r', tpc_comment.c_str() );
196 linadd( tnu->induc_dn * tnu->E2nu * EN1RYD * (*tnu->Pop),
197 22, chLabel.c_str() ,
'i', tpe_comment.c_str() );
212 multiplet_sum( multiplet, av_wl, sumxInt, sumxObsInt, sumcool, sumheat );
213 multiplet.resize( 0 );
216 " total emission in He-like intercombination lines from 2p3P to ground ");
230 for(
long ipLo=0; ipLo <
ipHe2p3P0; ipLo++ )
232 vector<long> EnterTheseLast;
233 for(
long ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
243 sp->
st[ipLo].l())==1 )
244 && (sp->
st[ipLo].l()>1)
245 && (sp->
st[ipHi].l()>1)
246 && ( sp->
st[ipHi].n() ==
250 if( (sp->
st[ipHi].S()==1)
251 && (sp->
st[ipLo].S()==1) )
255 else if( (sp->
st[ipHi].S()==3)
256 && (sp->
st[ipLo].S()==3) )
258 string comment_trans =
"";
265 HeTrList.push_back( sp->
trans(ipHi , ipLo+1) );
266 HeTrList.push_back( sp->
trans(ipHi+1, ipLo+1) );
268 if( multiplet.size() == 0 )
continue;
270 multiplet_sum( multiplet, av_wl, sumxInt, sumxObsInt, sumcool, sumheat );
271 multiplet.resize( 0 );
275 comment_trans =
"singlet to singlet: ";
277 comment_trans +=
"; ";
281 sp->
trans(ipHi+1,ipLo+1).
WLAng(), sumxInt, sumxObsInt, sumcool, sumheat,
282 comment_trans.c_str() );
289 HeTrList.push_back( sp->
trans(ipHi , ipLo) );
290 HeTrList.push_back( sp->
trans(ipHi+1, ipLo) );
292 if( multiplet.size() == 0 )
continue;
294 multiplet_sum( multiplet, av_wl, sumxInt, sumxObsInt, sumcool, sumheat );
295 multiplet.resize( 0 );
299 comment_trans =
"triplet to triplet: ";
301 comment_trans +=
"; ";
305 sp->
trans(ipHi,ipLo).
WLAng(), sumxInt, sumxObsInt, sumcool, sumheat,
306 comment_trans.c_str() );
310 else if( ipLo==
ipHe2s3S && ipHi == ipHe2p3P0 )
317 for(
long i=ipHe2p3P0; i <=
ipHe2p3P2; i++ )
333 HeTrList.push_back( sp->
trans(i, ipLo) );
336 if( multiplet.size() == 0 )
continue;
337 multiplet_sum( multiplet, av_wl, sumxInt, sumxObsInt, sumcool, sumheat );
338 multiplet.resize( 0 );
341 linadd(sumxObsInt, av_wl,
"TOTL",
'i',
342 "total emission in He-like lines, use average of three line wavelengths " );
347 string comment_trans =
"";
352 PutLineSum( sp->
trans(ipHi,ipLo), av_wl, sumxInt, sumxObsInt, sumcool, sumheat,
353 comment_trans.c_str() );
357 string comment_trans =
"";
365 linadd(sumxObsInt, av_wl, chLabel.c_str(),
'i',
366 "total emission in He-like lines, use average of three line wavelengths " );
379 string comment_trans =
"";
396 if( abs(
L_(ipHi) -
L_(ipLo) ) != 1 )
398 EnterTheseLast.push_back( ipHi );
405 string comment_trans =
"";
415 enum {DEBUG_LOC=
false};
418 if( nelem==1 && ipLo==0 && ipHi==1 )
431 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
433 string comment_trans =
"";
443 for(
long ipHi=
ipHe2p3P2+1; ipHi < nLoop; ipHi++ )
445 for(
long ipLo=ipHe2p3P0; ipLo <=
ipHe2p3P2; ++ipLo )
459 HeTrList.push_back( sp->
trans(ipHi, ipLo) );
464 multiplet.resize( 0 );
468 if( multiplet.size() == 0 )
continue;
469 multiplet_sum( multiplet, av_wl, sumxInt, sumxObsInt, sumcool, sumheat );
470 multiplet.resize( 0 );
472 string comment_trans =
"";
478 comment_trans.c_str() );
480 for(
long ipLo=
ipHe2p3P2+1; ipLo < nLoop-1; ipLo++ )
482 vector<long> EnterTheseLast;
483 for(
long ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
498 EnterTheseLast.push_back( ipHi );
502 string comment_trans =
"";
511 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
513 string comment_trans =
"";
531 for(
long ipDens = 0; ipDens <
NUMDENS; ++ipDens )
539 linadd( intens, wvlg, CaBLines[nelem][i].label,
'i',
"Case B intensity " );
555 double photons_3889_plus_7065 =
575 photons_3889_plus_7065 +=
582 linadd( photons_3889_plus_7065, 3889.,
"Pho+",
'i',
583 "photon sum given in Porter et al. 2007 (astro-ph/0611579)");
607 while (chLine[0] ==
'#');
617 # define chLine_LENGTH 1000
623 fprintf(
ioQQQ,
" lines_helium opening he1_case_b.dat\n");
625 ioDATA =
open_data(
"he1_case_b.dat",
"r" );
630 fprintf(
ioQQQ,
" lines_helium could not read first line of he1_case_b.dat.\n");
635 i1 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
636 i2 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
643 " lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
645 " lines_helium: I expected to find the number %i and got %li instead.\n" ,
647 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
652 if (
read_data_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
654 fprintf(
ioQQQ,
" lines_helium could not find line of densities in he1_case_ab.dat.\n");
660 for(
long j=0; !lgEOL && j <
NUMDENS; ++j)
666 if (
read_data_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
668 fprintf(
ioQQQ,
" lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
674 for (
long j=0; !lgEOL && j <
NUMTEMPS; ++j)
701 for(
long j = 0; j < i2; ++j )
705 CaBLines[nelem][j].
ipHi = -1;
706 CaBLines[nelem][j].
ipLo = -1;
707 strcpy( CaBLines[nelem][j].label ,
" " );
709 for(
long k = 0; k <
NUMDENS; ++k )
712 for(
long l = 0; l <
NUMTEMPS; ++l )
724 while(
read_data_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
727 if( (chLine[0] ==
' ') || (chLine[0]==
'\n') )
735 long ipLo = (long)
FFmtRead(chLine,&j,
sizeof(chLine),&lgEOL);
736 long ipHi = (long)
FFmtRead(chLine,&j,
sizeof(chLine),&lgEOL);
737 CaBLines[nelem][line].
ipLo = ipLo;
738 CaBLines[nelem][line].
ipHi = ipHi;
740 ASSERT( CaBLines[nelem][line].ipLo < CaBLines[nelem][line].ipHi );
742 strncpy( CaBLines[nelem][line].label,
"Ca B" , 4 );
743 CaBLines[nelem][line].
label[4] = 0;
751 for(
long ipDens = 0; ipDens <
NUMDENS; ++ipDens )
755 fprintf(
ioQQQ,
" lines_helium could not scan case B lines, current indices: %li %li\n",
756 CaBLines[nelem][line].ipHi,
757 CaBLines[nelem][line].ipLo );
761 char *chTemp = chLine;
763 long den = (long)
FFmtRead(chTemp,&j,
sizeof(chTemp),&lgEOL);
764 if( den != ipDens + 1 )
766 for(
long ipTe = 0; ipTe <
NUMTEMPS; ++ipTe )
769 if( (chTemp =
strchr_s( chTemp,
'\t' )) == NULL )
771 fprintf(
ioQQQ,
" lines_helium could not scan case B lines, current indices: %li %li\n",
772 CaBLines[nelem][line].ipHi,
773 CaBLines[nelem][line].ipLo );
777 sscanf( chTemp,
"%le" , &b );
793 STATIC double TempInterp2(
double* TempArray ,
double* ValueArray,
long NumElements,
double Te )
802 if( Te <= TempArray[0] )
804 return ValueArray[0] + log10( TempArray[0]/Te );
806 else if( Te >= TempArray[NumElements-1] )
808 return ValueArray[NumElements-1];
814 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
818 i0 =
max(
min(ipTe-ORDER/2,NumElements-ORDER-1),0);
819 rate =
lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, Te );
837 double dr_rate =
iso_sp[ipISO][nelem].
fb[i].DielecRecomb;
845 PutLine( tr,
"iso satellite line" );
static double **** CaBIntensity
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
string chIonLbl(const TransitionProxy &t)
NORETURN void TotalInsanity(void)
multi_arr< int, 3 > ipSatelliteLines
STATIC void randomize_inten(t_iso_sp *sp, long ipLo, long ipHi)
void set_xIntensity(const TransitionProxy &t)
string iso_comment_tran_levels(long ipISO, long nelem, long ipLo, long ipHi)
static double CaBDensities[NUMDENS]
double phots(const TransitionProxy &t)
long hunt_bisect(const T x[], long n, T xval)
t_iso_sp iso_sp[NISO][LIMELM]
STATIC void insert_trans(vector< TransitionProxy > &trList, TransitionProxy tr)
double xIonDense[LIMELM][LIMELM+1]
char * read_data_line(char *chLine, int size, FILE *ioDATA)
vector< two_photon > TwoNu
long int n_HighestResolved_max
realnum & EnergyWN() const
STATIC realnum get_multiplet_wl(vector< TransitionProxy > &multiHe, long ipHi, long ipLo)
STATIC void multiplet_sum(vector< TransitionProxy > &trList, realnum &av_wl, double &sum_inten, double &sum_obs_inten, double &sum_cool, double &sum_heat)
double & xIntensity() const
EmissionList::reference Emis() const
LinSv * linadd(double xEmiss, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
const char * strchr_s(const char *s, int c)
vector< vector< TransitionList > > SatelliteLines
STATIC void PutLineSum(const TransitionProxy &tr, const realnum av_wl, const double sumxInt, const double sumxObsInt, const double sumcool, const double sumheat, const char *chComment)
STATIC void GetStandardHeLines(void)
multi_arr< long, 3 > QuantumNumbers2Index
TransitionProxy trans(const long ipHi, const long ipLo)
double lagrange(const double x[], const double y[], long n, double xval)
STATIC void DoSatelliteLines(long nelem)
multi_arr< extra_tr, 2 > ex
void PutLine(const TransitionProxy &t, const char *chComment, const char *chLabelTemp, const ExtraInten &extra)
static double CaBTemps[NUMTEMPS]
CollisionProxy Coll() const
#define DEBUG_ENTRY(funcname)
double linint(const double x[], const double y[], long n, double xval)
int fprintf(const Output &stream, const char *format,...)
double & xObsIntensity() const
static stdLines ** CaBLines
vector< realnum > CaseBWlHeI
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
STATIC double TempInterp2(double *TempArray, double *ValueArray, long NumElements, double Te)
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
long int StuffComment(const char *chComment)