20 static void AulTTError(
const char chFilename[],
const char chLine[],
const char TT[] )
24 fprintf(
ioQQQ,
" PROBLEM File %s contains an invalid line.\n",chFilename);
25 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
26 fprintf(
ioQQQ,
" This transition already has an Aul value set for %s.\n",TT);
39 typename T::const_iterator i;
40 for( i = input.begin(); i != input.end(); i++ )
44 returnLev = i-input.begin();
58 for( DoubleLongPairVector::const_iterator i = dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
60 if( oldLev == i->second )
62 returnLev = i-dBaseStatesEnergy.begin();
84 static const int MAX_NUM_LEVELS = 999;
89 strcpy( chNRGFilename , chPrefix );
90 strcpy( chTPFilename , chPrefix );
91 strcpy( chCOLLFilename , chPrefix );
94 strcat( chNRGFilename ,
".nrg");
99 fprintf(
ioQQQ,
" atmdat_STOUT_readin opening %s:",chNRGFilename);
102 ioDATA =
open_data( chNRGFilename,
"r" );
111 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of %s.\n", chNRGFilename );
115 const long int iyr = 11, imo=10 , idy = 14;
116 long iyrread, imoread , idyread;
117 iyrread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
118 imoread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
119 idyread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
121 if(( iyrread != iyr ) ||
122 ( imoread != imo ) ||
126 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chNRGFilename );
128 " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
129 iyr, imo , idy ,iyrread, imoread , idyread );
135 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
138 if( chLine[0] !=
'#' && chLine[0] !=
'\n' && chLine[0] !=
'*' )
141 long n = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
146 else if( (chLine[0] ==
'*' && chLine[1] ==
'*' ) )
154 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
156 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not rewind %s.\n", chNRGFilename );
162 long HighestIndexInFile = nMolLevs;
184 long numMasterlist =
MIN2(
dBaseSpecies[intNS].numLevels_masterlist , HighestIndexInFile );
185 nMolLevs =
MAX2(nMolLevs,numMasterlist);
193 fprintf(
ioQQQ,
"Using STOUT spectrum %s (species: %s) with %li requested, only %li energy levels available.\n",
195 nMolLevs = HighestIndexInFile;
207 fprintf(
ioQQQ,
"Using STOUT spectrum %s (species: %s) with %li levels of %li available.\n",
208 dBaseSpecies[intNS].chLabel, chLabelChemical, nMolLevs, HighestIndexInFile );
221 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
232 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
234 for(
long ipLo = 0; ipLo < ipHi; ipLo++)
248 vector<double> dBaseStatesStwt(HighestIndexInFile,-1.0);
249 for(
long ii = 0; ii < HighestIndexInFile; ii++ )
251 dBaseStatesEnergy.push_back(make_pair(-1.0,ii));
258 fprintf(
ioQQQ,
"Number of Energy Levels in File: %li\n",HighestIndexInFile);
259 fprintf(
ioQQQ,
"Number of Energy Levels Cloudy is Currently Using: %li\n",nMolLevs);
260 fprintf(
ioQQQ,
"Species|File Index|Energy(WN)|StatWT\n");
264 bool lgSentinelReached =
false;
269 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of the energy level file.\n");
278 if( chLine[0] ==
'*' )
280 lgSentinelReached =
true;
284 if( chLine[0] !=
'#' )
287 if( chLine[0] ==
'\n')
290 index = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL) -1 ;
298 fprintf(
ioQQQ,
" PROBLEM An energy level index was less than 1 in file %s\n",chNRGFilename);
299 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
303 nrg = (double)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
304 stwt = (double)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
312 fprintf(
ioQQQ,
" PROBLEM End of line reached prematurely in file %s\n",chNRGFilename);
313 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
318 if( dBaseStatesEnergy.at(index).first == -1. && dBaseStatesStwt.at(index) == -1 )
320 dBaseStatesEnergy.at(index) = make_pair(nrg,index);
321 dBaseStatesStwt.at(index) = stwt;
325 fprintf(
ioQQQ,
" PROBLEM Duplicate energy level index in file %s\n",chNRGFilename);
326 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
331 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL );
332 if( !lgSentinelReached )
334 fprintf(
ioQQQ,
" PROBLEM End of data sentinel was not reached in file %s\n",chNRGFilename);
335 fprintf(
ioQQQ,
" Check that stars (*****) appear after the last line of data and start in the first column of that line.\n");
341 sort(dBaseStatesEnergy.begin(),dBaseStatesEnergy.end());
345 fprintf(
ioQQQ,
"\n\n***Energy levels have been sorted in order of increasing energy.***\n");
346 fprintf(
ioQQQ,
"Species|File Index|Sorted Index|Energy(WN)|StatWT\n");
350 for(DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin();i!=dBaseStatesEnergy.end();i++)
352 long oldindex = i->second;
353 long index = i - dBaseStatesEnergy.begin();
354 double nrg = i->first;
355 double stwt = dBaseStatesStwt.at(oldindex);
357 if( index >= nMolLevs )
365 dBaseStates[intNS][index].energy().set(nrg,
"cm^-1");
371 double fenergyWN = 0.;
375 int ipHi = (*tr).ipHi();
376 int ipLo = (*tr).ipLo();
381 long ipLoFile = dBaseStatesEnergy[ipLo].second;
382 long ipHiFile = dBaseStatesEnergy[ipHi].second;
387 if( fenergyWN == 0. &&
dBaseStates[intNS][ipHi].g() !=
dBaseStates[intNS][ipLo].g() && abs(ipHiFile - ipLoFile) == 1)
393 fprintf(
ioQQQ,
"\nCaution: A %s transition between adjacent energy levels has zero energy.\n",
dBaseSpecies[intNS].chLabel);
394 fprintf(
ioQQQ,
"The levels may have been sorted by energy since being read from the data files.\n");
395 fprintf(
ioQQQ,
"The sorted lower and upper levels are %i and %i.\n",ipLo+1,ipHi+1);
396 fprintf(
ioQQQ,
"The level indices as they appear in the Stout energy level data file, %s, are %li and %li.\n",chNRGFilename,ipLoFile+1,ipHiFile+1);
397 fprintf(
ioQQQ,
"Verify with the data source that the correct energies are listed in the Stout data file.\n");
403 fprintf(
ioQQQ,
"The levels may have been sorted by energy since being read from the data files.\n");
404 fprintf(
ioQQQ,
"The sorted lower and upper levels are %i and %i.\n",ipLo+1,ipHi+1);
405 fprintf(
ioQQQ,
"The level indices as they appear in the Stout energy level data file, %s, are %li and %li.\n",chNRGFilename,ipLoFile+1,ipHiFile+1);
406 fprintf(
ioQQQ,
"Check the data file for missing or duplicate energies.\n");
412 (*tr).WLAng() = (
realnum)(1e+8/(*tr).EnergyWN()/
RefIndex((*tr).EnergyWN()));
419 strcat( chTPFilename ,
".tp");
427 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of %s.\n", chTPFilename );
431 iyrread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
432 imoread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
433 idyread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
435 if(( iyrread != iyr ) ||
436 ( imoread != imo ) ||
440 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chTPFilename );
442 " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
443 iyr, imo , idy ,iyrread, imoread , idyread );
449 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
452 if( chLine[0] !=
'#' && chLine[0] !=
'\n' && chLine[0] !=
'*' )
455 long n = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
460 else if( (chLine[0] ==
'*' && chLine[1] ==
'*' ) )
467 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
469 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not rewind %s.\n", chTPFilename );
479 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of the transition probability file.\n");
484 lgSentinelReached =
false;
485 static const int intNumCols = 6;
500 bool ***lgLineStrengthTT;
501 lgLineStrengthTT = (
bool ***)
MALLOC(nMolLevs *
sizeof(
bool**));
502 for(
int ii = 0; ii < nMolLevs; ii++ )
504 lgLineStrengthTT[ii] = (
bool **)
MALLOC(nMolLevs *
sizeof(
bool*));
505 for(
int j = 0; j < nMolLevs; j++ )
507 lgLineStrengthTT[ii][j] = (
bool *)
MALLOC(intNumCols *
sizeof(
bool));
512 for(
int ii = 0; ii < nMolLevs; ii++ )
514 for(
int j = 0; j < nMolLevs; j++ )
516 for(
int k = 0; k < intNumCols; k++ )
518 lgLineStrengthTT[ii][j][k] =
false;
526 fprintf(
ioQQQ,
"Radiative Data File: %s\n",chTPFilename);
527 fprintf(
ioQQQ,
"Species|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Data Type (A,G,S)|Data\n");
533 if( chLine[0] ==
'*' )
535 lgSentinelReached =
true;
540 if( chLine[0] !=
'#' )
543 if( chLine[0] ==
'\n')
553 long ipLoInFile = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
554 long ipHiInFile = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
555 long ipLo = ipLoInFile - 1;
556 long ipHi = ipHiInFile - 1;
557 tpData = (double)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
568 if( ipLo < 0 || ipLo >= nMolLevs || ipHi < 0 || ipHi >= nMolLevs )
588 ipLo+1,ipHi+1,chLine[0],tpData);
593 fprintf(
ioQQQ,
" PROBLEM End of line reached prematurely in file %s\n",chTPFilename);
594 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
602 if( !(*tr).hasEmis() )
604 (*tr).AddLine2Stack();
605 (*tr).Emis().Aul() = 0.;
606 (*tr).Emis().gf() = 0.;
614 (*tr).Emis().Aul() += tpData;
616 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
620 else if(
nMatch(
"G",chLine) )
624 (*tr).Emis().gf() += tpData;
626 (*tr).Emis().Aul() = (
realnum)
eina((*tr).Emis().gf(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
629 else if(
nMatch(
"S",chLine) )
631 static const double BOHR_MAGNETON = ELEM_CHARGE_ESU*H_BAR/2/ELECTRON_MASS/SPEEDLIGHT;
637 if( lgLineStrengthTT[ipLo][ipHi][0] )
644 static const double E1Coeff = 64*
powi(PI,4)*
pow(ELEM_CHARGE_ESU,2)*
pow2(BOHR_RADIUS_CM)/3/HPLANCK/
pow3(1e-8);
647 (*tr).Emis().Aul() += E1Coeff*tpData/(*(*tr).Hi()).g()/
pow3((*tr).EnergyAng());
648 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
649 lgLineStrengthTT[ipLo][ipHi][0]=
true;
651 else if(
nMatch(
"E2",chLine) )
653 if( lgLineStrengthTT[ipLo][ipHi][1] )
660 static const double E2Coeff = 64*
powi(PI,6)*
pow2(ELEM_CHARGE_ESU)*
pow4(BOHR_RADIUS_CM)/15/HPLANCK/
powi(1e-8,5);
663 (*tr).Emis().Aul() += E2Coeff*tpData/(*(*tr).Hi()).g()/
powi((*tr).EnergyAng(),5);
664 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
665 lgLineStrengthTT[ipLo][ipHi][1]=
true;
667 else if(
nMatch(
"E3",chLine) )
669 if( lgLineStrengthTT[ipLo][ipHi][2] )
676 static const double E3Coeff = 2048*
powi(PI,8)*
pow2(ELEM_CHARGE_ESU)*
powi(BOHR_RADIUS_CM,6)/4725/HPLANCK/
powi(1e-8,7);
681 (*tr).Emis().Aul() += E3Coeff*tpData/(*(*tr).Hi()).g()/
powi((*tr).EnergyAng(),7);
682 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
683 lgLineStrengthTT[ipLo][ipHi][2]=
true;
685 else if(
nMatch(
"M1",chLine) )
687 if( lgLineStrengthTT[ipLo][ipHi][3] )
694 static const double M1Coeff = 64*
powi(PI,4)*
pow2(BOHR_MAGNETON)/3/HPLANCK/
pow3(1e-8);
697 (*tr).Emis().Aul() += M1Coeff*tpData/(*(*tr).Hi()).g()/
pow3((*tr).EnergyAng());
698 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
699 lgLineStrengthTT[ipLo][ipHi][3]=
true;
701 else if(
nMatch(
"M2",chLine) )
703 if( lgLineStrengthTT[ipLo][ipHi][4] )
710 static const double M2Coeff = 64*
powi(PI,6)*
pow2(BOHR_MAGNETON)*
pow2(BOHR_RADIUS_CM)/15/HPLANCK/
powi(1e-8,5);
715 (*tr).Emis().Aul() += M2Coeff*tpData/(*(*tr).Hi()).g()/
powi((*tr).EnergyAng(),5);
716 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
717 lgLineStrengthTT[ipLo][ipHi][4]=
true;
719 else if(
nMatch(
"M3",chLine) )
721 if( lgLineStrengthTT[ipLo][ipHi][5] )
728 static const double M3Coeff = 2048*
powi(PI,8)*
pow2(BOHR_MAGNETON)*
pow4(BOHR_RADIUS_CM)/4725/HPLANCK/
powi(1e-8,7);
733 (*tr).Emis().Aul() += M3Coeff*tpData/(*(*tr).Hi()).g()/
powi((*tr).EnergyAng(),7);
734 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
735 lgLineStrengthTT[ipLo][ipHi][5]=
true;
739 fprintf(
ioQQQ,
" PROBLEM File %s contains an invalid line.\n",chTPFilename);
740 fprintf(
ioQQQ,
" The line strength does not list a valid transition type.\n");
741 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
756 fprintf(
ioQQQ,
" PROBLEM File %s contains an invalid line.\n",chTPFilename);
757 fprintf(
ioQQQ,
" Data must either be in the form of Aul, gf, or S(line strength) and list"
758 " the data type in the first colum as A, G, or S respectively.");
759 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
764 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL );
765 if( !lgSentinelReached )
767 fprintf(
ioQQQ,
" PROBLEM End of data sentinel was not reached in file %s\n",chTPFilename);
768 fprintf(
ioQQQ,
" Check that stars (*****) appear after the last line of data and start in the first column of that line.");
777 strcat( chCOLLFilename ,
".coll");
780 ioDATA =
open_data( chCOLLFilename,
"r" );
785 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of %s.\n", chCOLLFilename );
789 iyrread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
790 imoread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
791 idyread = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
793 if(( iyrread != iyr ) ||
794 ( imoread != imo ) ||
798 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chCOLLFilename );
800 " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
801 iyr, imo , idy ,iyrread, imoread , idyread );
811 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of the collision data file.\n");
817 for(
long ipHi=0; ipHi<nMolLevs; ipHi++ )
819 for(
long ipLo=0; ipLo<nMolLevs; ipLo++ )
830 vector<double> temps;
831 long ipCollider = -1;
832 lgSentinelReached =
false;
837 fprintf(
ioQQQ,
"Collision Data File: %s\n",chCOLLFilename);
839 fprintf(
ioQQQ,
"Species|Data Type (CS,RATE)|Collider|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Data\n");
846 if( chLine[0] ==
'*' )
848 lgSentinelReached =
true;
853 if( chLine[0] !=
'#' )
858 if( chLine[0] ==
'\n')
865 if(
nMatch(
"TEMP",chLine) )
878 FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
885 temps.resize(numpoints);
886 for(
int j = 0; j < numpoints; j++ )
892 for(
int j = 0; j < numpoints; j++ )
894 temps[j] = (double)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
911 if(
nMatch(
"RATE", chLine) )
914 if(
nMatch(
"ELECTRON",chLine ) )
918 else if(
nMatch(
"PROTON",chLine ) ||
nMatch(
"H+",chLine ) )
922 else if(
nMatch(
"HE+2", chLine ) )
926 (void)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
928 else if(
nMatch(
"HE+ ",chLine ) ||
nMatch(
"HE+\t", chLine) )
932 else if(
nMatch(
"H2 ",chLine ) ||
nMatch(
"H2\t",chLine ) )
935 (void)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
937 if(
nMatch(
"ORTHO",chLine ) )
941 else if(
nMatch(
"PARA",chLine ) )
950 else if(
nMatch(
"HE ",chLine ) ||
nMatch(
"HE\t",chLine ) )
954 else if(
nMatch(
"H ",chLine ) ||
nMatch(
"H\t",chLine ) )
960 fprintf(
ioQQQ,
" PROBLEM atmdat_STOUT_readin: Each line of the collision data"
961 "file must specify the collider.\n");
962 fprintf(
ioQQQ,
" Possible colliders are: ELECTRON, PROTON, HE, H,"
963 " HE+2, HE+, H2, H2 ORTHO, H2 PARA\n");
969 fprintf(
ioQQQ,
" PROBLEM atmdat_STOUT_readin: The collision "
970 "file must specify temperatures before the collision strengths");
974 long ipLoInFile = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
975 long ipHiInFile = (long)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
976 long ipLo = ipLoInFile-1;
977 long ipHi = ipHiInFile-1;
983 if( ipLo < 0 || ipLo >= nMolLevs || ipHi < 0 || ipHi >= nMolLevs )
1002 dBaseSpecies[intNS].chLabel,isRate?
"RATE":
"CS",ipCollider,
1003 ipLoInFile,ipHiInFile,ipLo+1,ipHi+1);
1007 StoutCollData[intNS].lgIsRate(ipHi,ipLo,ipCollider) = isRate;
1010 StoutCollData[intNS].setpoints(ipHi,ipLo,ipCollider,numpoints);
1013 for(
int j = 0; j < numpoints; j++ )
1015 StoutCollData[intNS].temps(ipHi,ipLo,ipCollider)[j] = temps[j];
1016 StoutCollData[intNS].collstrs(ipHi,ipLo,ipCollider)[j] = (double)
FFmtRead(chLine,&ipFFmt,
sizeof(chLine),&lgEOL);
1021 if(
StoutCollData[intNS].temps(ipHi,ipLo,ipCollider)[j] <= 0 ||
1022 StoutCollData[intNS].collstrs(ipHi,ipLo,ipCollider)[j] <= 0 )
1024 fprintf(
ioQQQ,
"PROBLEM: A Stout temperature or collision strength is less than or equal to zero.\n");
1026 fprintf(
ioQQQ,
"numpoint = %i\tTemp = %e\tCS = %e\n",j,temps[j],
1036 fprintf(
ioQQQ,
" PROBLEM File %s contains an invalid line.\n",chCOLLFilename);
1037 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
1043 fprintf(
ioQQQ,
" PROBLEM End of line reached prematurely in file %s\n",chCOLLFilename);
1044 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
1049 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL );
1050 if( !lgSentinelReached )
1052 fprintf(
ioQQQ,
" PROBLEM End of data sentinel was not reached in file %s\n",chCOLLFilename);
1053 fprintf(
ioQQQ,
" Check that stars (*****) appear after the last line of data and start in the first column.");
1059 for(
int ii = 0; ii < nMolLevs; ii++ )
1061 for(
int j = 0; j < nMolLevs; j++ )
1063 free( lgLineStrengthTT[ii][j] );
1065 free( lgLineStrengthTT[ii] );
1067 free( lgLineStrengthTT );
1076 int intsplinepts,intTranType,intxs;
1077 long int nMolLevs,nMolExpLevs,nElvlcLines,nTheoLevs;
1078 FILE *ioElecCollData=NULL, *ioProtCollData=NULL;
1079 realnum fstatwt,fenergyWN,fWLAng,fenergy,feinsteina;
1080 double fScalingParam,fEnergyDiff;
1081 const char chCommentChianti =
'#';
1089 bool lgProtonData=
false;
1092 static const int MAX_NUM_LEVELS = 999;
1097 strcpy( chEnFilename , chPrefix );
1098 strcpy( chTraFilename , chPrefix );
1099 strcpy( chEleColFilename , chPrefix );
1100 strcpy( chProColFilename , chPrefix );
1104 strcat( chEnFilename ,
".elvlc");
1109 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chEnFilename);
1111 fstream elvlcstream;
1115 strcat( chTraFilename ,
".wgfa");
1120 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chTraFilename);
1126 strcat( chEleColFilename ,
".splups");
1127 uncaps( chEleColFilename );
1131 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chEleColFilename);
1133 ioElecCollData =
open_data( chEleColFilename,
"r" );
1136 strcat( chProColFilename ,
".psplups");
1137 uncaps( chProColFilename );
1141 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chProColFilename);
1146 lgProtonData =
true;
1150 lgProtonData =
false;
1155 const int eof_col = 5;
1157 const int lvl_nrg_col=16;
1159 const int lvl_skipto_nrg = 40;
1161 const int lvl_eof_to_nrg = lvl_skipto_nrg - eof_col + 1;
1163 const int lvl_skip_ryd = 15;
1167 if (elvlcstream.is_open())
1170 char otemp[eof_col];
1171 char exptemp[lvl_nrg_col],theotemp[lvl_nrg_col];
1172 double tempexpenergy = 0.,theoenergy = 0.;
1177 elvlcstream.get(otemp,eof_col);
1183 elvlcstream.seekg(lvl_eof_to_nrg,ios::cur);
1184 elvlcstream.get(exptemp,lvl_nrg_col);
1185 tempexpenergy = (double) atof(exptemp);
1186 if( tempexpenergy != 0.)
1189 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
1190 elvlcstream.get(theotemp,lvl_nrg_col);
1191 theoenergy = (double) atof(theotemp);
1192 if( theoenergy != 0. )
1195 elvlcstream.ignore(INT_MAX,
'\n');
1198 elvlcstream.seekg(0,ios::beg);
1203 bool lgChiaBadTheo =
false;
1206 lgChiaBadTheo =
true;
1209 fprintf(
ioQQQ,
"Switching to the experimental levels for this species.");
1212 long HighestIndexInFile = -1;
1217 HighestIndexInFile = nMolExpLevs;
1221 HighestIndexInFile = nElvlcLines;
1224 dBaseSpecies[intNS].numLevels_max = HighestIndexInFile;
1240 fprintf(
ioQQQ,
"The number of energy levels is non-positive in datafile %s.\n", chEnFilename );
1246 long numMasterlist =
MIN2(
dBaseSpecies[intNS].numLevels_masterlist , HighestIndexInFile );
1247 nMolLevs =
MAX2(nMolLevs,numMasterlist);
1251 if (
dBaseSpecies[intNS].setLevels > HighestIndexInFile)
1255 fprintf(
ioQQQ,
"Using CHIANTI spectrum %s (species: %s) with %li requested, only %li energy levels available.\n",
1257 nMolLevs = HighestIndexInFile;
1274 fprintf(
ioQQQ,
"Using CHIANTI spectrum %s (species: %s) with %li experimental energy levels of %li available.\n",
1275 dBaseSpecies[intNS].chLabel, chLabelChemical, nMolLevs , nMolExpLevs );
1281 fprintf(
ioQQQ,
"Using CHIANTI spectrum %s (species: %s) with %li theoretical energy levels of %li available.\n",
1282 dBaseSpecies[intNS].chLabel, chLabelChemical, nMolLevs , nElvlcLines );
1291 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
1300 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
1302 for(
long ipLo = 0; ipLo < ipHi; ipLo++)
1318 vector<long> intExperIndex(nElvlcLines,-1);
1321 vector<double> dBaseStatesStwt(HighestIndexInFile,-1.0);
1322 for(
long ii = 0; ii < HighestIndexInFile; ii++ )
1324 dBaseStatesEnergy.push_back(make_pair(-1.0,ii));
1328 const int lvl_skipto_statwt = 37;
1330 const int lvl_statwt_col = 4;
1334 for(
long ipLev=0; ipLev<nElvlcLines; ipLev++ )
1336 if(elvlcstream.is_open())
1338 char gtemp[lvl_statwt_col],thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
1339 elvlcstream.seekg(lvl_skipto_statwt,ios::cur);
1340 elvlcstream.get(gtemp,lvl_statwt_col);
1341 fstatwt = (
realnum)atof(gtemp);
1342 elvlcstream.get(thtemp,lvl_nrg_col);
1343 fenergy = (double) atof(thtemp);
1347 fprintf(
ioQQQ,
" WARNING: A positive non zero value is expected for the "
1348 "statistical weight but was not found in %s"
1349 " level %li\n", chEnFilename,ipLev);
1361 if( fenergy != 0. || ipLev == 0 )
1363 dBaseStatesEnergy.at(ncounter).first = fenergy;
1364 dBaseStatesEnergy.at(ncounter).second = ncounter;
1365 dBaseStatesStwt.at(ncounter) = fstatwt;
1366 intExperIndex.at(ipLev) = ncounter;
1371 intExperIndex.at(ipLev) = -1;
1376 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
1377 elvlcstream.get(obtemp,lvl_nrg_col);
1378 fenergy = (double) atof(obtemp);
1379 if(fenergy != 0. || ipLev == 0)
1381 dBaseStatesEnergy.at(ipLev).first = fenergy;
1382 dBaseStatesEnergy.at(ipLev).second = ipLev;
1383 dBaseStatesStwt.at(ipLev) = fstatwt;
1387 dBaseStatesEnergy.at(ipLev).first = -1.;
1388 dBaseStatesEnergy.at(ipLev).second = ipLev;
1389 dBaseStatesStwt.at(ipLev) = -1.;
1393 elvlcstream.ignore(INT_MAX,
'\n');
1397 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEnFilename);
1398 fclose( ioProtCollData );
1403 elvlcstream.close();
1409 for( vector<long>::iterator i = intExperIndex.begin(); i != intExperIndex.end(); i++ )
1412 long iPrt = (i-intExperIndex.begin())+1;
1416 for( DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
1419 long iPrt = (i-dBaseStatesEnergy.begin())+1;
1421 (i->second)+1,i->first,dBaseStatesStwt.at(i->second));
1426 sort(dBaseStatesEnergy.begin(),dBaseStatesEnergy.end());
1430 for( DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
1433 long iPrt = (i-dBaseStatesEnergy.begin())+1;
1434 if( iPrt > nMolLevs )
1437 (i->second)+1,i->first,dBaseStatesStwt.at(i->second));
1444 fprintf(
ioQQQ,
"Number of Experimental Energy Levels in File: %li\n",nMolExpLevs);
1448 fprintf(
ioQQQ,
"Number of Theoretical Energy Levels in File: %li\n",nElvlcLines);
1450 fprintf(
ioQQQ,
"Number of Energy Levels Cloudy is Currently Using: %li\n",nMolLevs);
1451 fprintf(
ioQQQ,
"Species|File Index|Cloudy Index|StatWT|Energy(WN)\n");
1454 for( DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
1457 long ipLevNew = i - dBaseStatesEnergy.begin();
1458 long ipLevFile = -1;
1460 if( ipLevNew >= nMolLevs )
1469 ipLevFile = i->second;
1477 dBaseStates[intNS][ipLevNew].g() = dBaseStatesStwt.at(i->second);
1478 dBaseStates[intNS][ipLevNew].energy().set(i->first,
"cm^-1");
1493 int ipHi = (*tr).ipHi();
1494 int ipLo = (*tr).ipLo();
1497 (*tr).EnergyWN() = fenergyWN;
1509 if (wgfastream.is_open())
1512 char otemp[eof_col];
1515 wgfastream.get(otemp,eof_col);
1516 wgfastream.ignore(INT_MAX,
'\n');
1517 if( otemp[0] == chCommentChianti )
continue;
1521 wgfastream.seekg(0,ios::beg);
1524 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chTraFilename);
1529 fprintf(
ioQQQ,
"\n\nTransition Probability File: %s\n",chTraFilename);
1530 fprintf(
ioQQQ,
"Species|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Wavelength(A)|Ein A\n");
1535 const int line_index_col = 6;
1537 const int line_nrg_to_eina =15;
1539 const int line_eina_col = 16;
1540 char lvltemp[line_index_col];
1542 if (wgfastream.is_open())
1544 for (
long ii = 0;ii<wgfarows;ii++)
1546 wgfastream.get(lvltemp,line_index_col);
1547 if( lvltemp[0] == chCommentChianti )
1549 wgfastream.ignore(INT_MAX,
'\n');
1553 long ipLoInFile = atoi(lvltemp);
1554 wgfastream.get(lvltemp,line_index_col);
1555 long ipHiInFile = atoi(lvltemp);
1558 long ipLo = ipLoInFile - 1;
1559 long ipHi = ipHiInFile - 1;
1566 if( intExperIndex[ipLo] == -1 || intExperIndex[ipHi] == -1 )
1568 wgfastream.ignore(INT_MAX,
'\n');
1579 long testlo = -1, testhi = -1;
1586 catch ( out_of_range& )
1590 fprintf(
ioQQQ,
"NOTE: An out of range exception has occurred"
1591 " reading in data from %s\n",chTraFilename);
1592 fprintf(
ioQQQ,
" The line in the file containing the unidentifiable"
1593 " levels has been ignored.\n");
1595 " This message is just for documentation.\n");
1598 wgfastream.ignore(INT_MAX,
'\n');
1602 if( testlo == -1 || testhi == -1 )
1604 wgfastream.ignore(INT_MAX,
'\n');
1614 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1617 wgfastream.ignore(INT_MAX,
'\n');
1623 fprintf(
ioQQQ,
" WARNING: Upper level = lower for a radiative transition in %s. Ignoring.\n", chTraFilename );
1624 wgfastream.ignore(INT_MAX,
'\n');
1649 char trantemp[lvl_nrg_col];
1650 wgfastream.get(trantemp,lvl_nrg_col);
1651 fWLAng = (
realnum)atof(trantemp);
1683 wgfastream.seekg(line_nrg_to_eina,ios::cur);
1684 wgfastream.get(trantemp,line_eina_col);
1685 feinsteina = (
realnum)atof(trantemp);
1686 if( feinsteina == 0. )
1688 static bool notPrintedYet =
true;
1691 fprintf(
ioQQQ,
" WARNING: Radiative rate(s) equal to zero in %s.\n", chTraFilename );
1692 notPrintedYet =
false;
1694 wgfastream.ignore(INT_MAX,
'\n');
1702 fixit(
"may need to do something with these rates");
1704 wgfastream.getline(chLine,INT_MAX);
1706 if(
nMatch(
"auto", chLine) )
1708 if( (*tr).hasEmis() )
1710 (*tr).Emis().AutoIonizFrac() =
1711 feinsteina/((*tr).Emis().Aul() + feinsteina);
1712 ASSERT( (*tr).Emis().AutoIonizFrac() >= 0. );
1713 ASSERT( (*tr).Emis().AutoIonizFrac() <= 1. );
1718 if( (*tr).hasEmis() )
1720 fprintf(
ioQQQ,
" PROBLEM duplicate transition found by atmdat_chianti in %s, "
1721 "wavelength=%f\n", chTraFilename,fWLAng);
1725 (*tr).AddLine2Stack();
1726 (*tr).Emis().Aul() = feinsteina;
1728 fenergyWN = (
realnum)(1e+8/fWLAng);
1732 (*tr).EnergyWN() = fenergyWN;
1734 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
1739 else fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chTraFilename);
1744 for(
long ipHi=0; ipHi<nMolLevs; ipHi++ )
1747 for(
long ipLo=0; ipLo<nMolLevs; ipLo++ )
1754 for(
long ipHi=0; ipHi<nMolLevs; ipHi++ )
1756 for(
long ipLo=0; ipLo<nMolLevs; ipLo++ )
1776 for(
long ipCollider=0; ipCollider<=1; ipCollider++ )
1782 strcpy( chFilename, chEleColFilename );
1788 strcpy( chFilename, chProColFilename );
1798 fprintf(
ioQQQ,
"\n\nCollision Data File: %s\n",chTraFilename);
1799 fprintf(
ioQQQ,
"Species|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Spline Points\n");
1802 fstream splupsstream;
1806 const int cs_eof_col = 4;
1808 const int cs_index_col = 4;
1810 const int cs_trantype_col = 4;
1813 const int cs_values_col = 11;
1815 if (splupsstream.is_open())
1819 long splupslines = -1;
1820 char otemp[cs_eof_col];
1823 splupsstream.get(otemp,cs_eof_col);
1824 splupsstream.ignore(INT_MAX,
'\n');
1828 splupsstream.seekg(0,ios::beg);
1830 for (
int m = 0;m<splupslines;m++)
1834 splupsstream.seekg(6,ios::cur);
1838 splupsstream.get(otemp,cs_index_col);
1839 long ipLo = atoi(otemp)-1;
1840 splupsstream.get(otemp,cs_index_col);
1841 long ipHi = atoi(otemp)-1;
1843 long ipLoFile = ipLo;
1844 long ipHiFile = ipHi;
1851 if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1853 splupsstream.ignore(INT_MAX,
'\n');
1864 long testlo = -1, testhi = -1;
1874 catch ( out_of_range& )
1878 fprintf(
ioQQQ,
"NOTE: An out of range exception has occurred"
1879 " reading in data from %s\n",chEleColFilename);
1880 fprintf(
ioQQQ,
" The line in the file containing the unidentifiable"
1881 " levels has been ignored.\n");
1883 " This message is for documentation.\n");
1885 splupsstream.ignore(INT_MAX,
'\n');
1889 if( testlo == -1 || testhi == -1 )
1891 splupsstream.ignore(INT_MAX,
'\n');
1901 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1904 splupsstream.ignore(INT_MAX,
'\n');
1921 splupsstream.get(otemp,cs_trantype_col);
1922 intTranType = atoi(otemp);
1923 char qtemp[cs_values_col];
1924 splupsstream.get(qtemp,cs_values_col);
1926 splupsstream.get(qtemp,cs_values_col);
1927 fEnergyDiff = atof(qtemp);
1929 splupsstream.get(qtemp,cs_values_col);
1930 fScalingParam = atof(qtemp);
1933 ASSERT( ipLo >= 0 && ipLo < nMolLevs );
1934 ASSERT( ipHi >= 0 && ipHi < nMolLevs );
1938 const int CHIANTI_SPLINE_MAX=9, CHIANTI_SPLINE_MIN=5;
1943 (
double *)
MALLOC((
unsigned long)(CHIANTI_SPLINE_MAX)*
sizeof(
double));
1945 (
double *)
MALLOC((
unsigned long)(CHIANTI_SPLINE_MAX)*
sizeof(
double));
1948 for( intsplinepts=0; intsplinepts<=CHIANTI_SPLINE_MAX; intsplinepts++ )
1951 char p = splupsstream.peek();
1958 if( intsplinepts >= CHIANTI_SPLINE_MAX )
1960 fprintf(
ioQQQ,
" WARNING: More spline points than expected in %s, indices %3li %3li. Ignoring extras.\n", chFilename, ipHi, ipLo );
1963 ASSERT( intsplinepts < CHIANTI_SPLINE_MAX );
1966 splupsstream.get(qtemp,cs_values_col);
1974 if( intTranType < 6 )
1975 temp =
max( temp, 0. );
1976 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intsplinepts] = temp;
1985 ASSERT( intsplinepts > 2 );
1994 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].ScalingParam = fScalingParam;
1998 vector<double> xs (intsplinepts),
2002 for(intxs=0;intxs<intsplinepts;intxs++)
2004 double coeff = (double)1/(intsplinepts-1);
2005 xs[intxs] = coeff*intxs;
2006 spl[intxs] =
AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
2009 spline(&xs[0], &spl[0],intsplinepts,2e31,2e31,&spl2[0]);
2012 for(
long ii=0; ii<intsplinepts; ii++)
2014 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer[ii] = spl2[ii];
2017 splupsstream.ignore(INT_MAX,
'\n');
2019 splupsstream.close();
2024 fclose( ioElecCollData );
2026 fclose( ioProtCollData );
long getSortedLevel(const T &input, long oldLev)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
vector< StoutCollArray > StoutCollData
static const bool DEBUGSTATE
const int FILENAME_PATH_LENGTH_2
NORETURN void TotalInsanity(void)
double eina(double gf, double enercm, double gup)
long nMatch(const char *chKey, const char *chCard)
vector< pair< double, long > > DoubleLongPairVector
vector< multi_arr< int, 2 > > ipdBaseTrans
void atmdat_STOUT_readin(long intNS, char *chFileName)
double RefIndex(double EnergyWN)
STATIC void spectral_to_chemical(char *chLabelChemical, char *chLabel, long &nelem, long &IonStg)
void uncaps(char *chCard)
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
const ios_base::openmode mode_r
static const double aulThreshold
void swap(count_ptr< T > &a, count_ptr< T > &b)
double energy(const genericState &gs)
static void AulTTError(const char chFilename[], const char chLine[], const char TT[])
void setProperties(species &sp)
double powi(double, long int)
vector< multi_arr< CollSplinesArray, 3 > > AtmolCollSplines
const long nDefaultPhotoLevelsFe
double GetGF(double trans_prob, double enercm, double gup)
string db_comment_tran_levels(long ipLoFile, long ipHiFile)
#define DEBUG_ENTRY(funcname)
vector< qList > dBaseStates
const double ENERGY_MIN_WN
vector< species > dBaseSpecies
int fprintf(const Output &stream, const char *format,...)
void atmdat_CHIANTI_readin(long intNS, char *chFileName)
double pow(double x, int i)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
vector< TransitionList > dBaseTrans
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)