27 static const int MAX_NUM_LEVELS = 999;
33 strcpy( chNRGFilename , chPrefix );
34 strcpy( chTPFilename , chPrefix );
35 strcpy( chCOLLFilename , chPrefix );
38 strcat( chNRGFilename ,
".nrg");
41 printf(
"The data file is %s \n",chNRGFilename);
45 fprintf(
ioQQQ,
" atmdat_STOUT_readin opening %s:",chNRGFilename);
58 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of %s.\n", chNRGFilename );
62 const long int iyr = 11, imo=10 , idy = 14;
63 long iyrread, imoread , idyread;
64 iyrread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
65 imoread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
66 idyread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
68 if(( iyrread != iyr ) ||
73 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chNRGFilename );
75 " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
76 iyr, imo , idy ,iyrread, imoread , idyread );
82 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
85 if( chLine[0] !=
'#' && chLine[0] !=
'\n' && chLine[0] !=
'*' )
88 long n = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
96 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
98 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not rewind %s.\n", chNRGFilename );
104 long HighestIndexInFile = nMolLevs;
106 nMolLevs =
MIN3(
atmdat.nStoutMaxLevels,nMolLevs, MAX_NUM_LEVELS );
108 if(
atmdat.lgStoutPrint ==
true)
110 fprintf(
ioQQQ,
" Using STOUT model %s with %li levels of %li available.\n",
111 dBaseSpecies[intNS].chLabel , nMolLevs , HighestIndexInFile );
124 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
134 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
136 for(
long ipLo = 0; ipLo < ipHi; ipLo++)
148 bool lgSentinelReached =
false;
153 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of the energy level file.\n");
162 if( chLine[0] ==
'*' )
164 lgSentinelReached =
true;
168 if( chLine[0] !=
'#' )
171 if( chLine[0] ==
'\n')
174 index = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL) -1 ;
178 fprintf(
ioQQQ,
" PROBLEM An energy level index was less than 1 in file %s\n",chNRGFilename);
179 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
183 nrg = (double)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
184 stwt = (double)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
188 fprintf(
ioQQQ,
" PROBLEM End of line reached prematurely in file %s\n",chNRGFilename);
189 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
193 if( index >= nMolLevs )
199 dBaseStates[intNS][index].energy().set(nrg,
"cm^-1");
208 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL );
209 if( !lgSentinelReached )
211 fprintf(
ioQQQ,
" PROBLEM End of data sentinel was not reached in file %s\n",chNRGFilename);
212 fprintf(
ioQQQ,
" Check that stars (*****) appear after the last line of data and start in the first column of that line.\n");
220 printf(
"DEBUG STOUT ENERGY LEVELS:\n");
221 for(
int i = 0; i< nMolLevs; i++ )
228 double fenergyWN = 0.;
232 int ipHi = (*tr).ipHi();
233 int ipLo = (*tr).ipLo();
238 printf(
"\nWARNING: The %s transition between level %i and %i has zero or negative energy.\n",
240 printf(
"Check the Stout energy level data file (*.nrg) of this species for missing or duplicate energies.\n");
244 (*tr).WLAng() = (
realnum)(1e+8/(*tr).EnergyWN()/
RefIndex((*tr).EnergyWN()));
251 strcat( chTPFilename ,
".tp");
259 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of %s.\n", chTPFilename );
263 iyrread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
264 imoread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
265 idyread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
267 if(( iyrread != iyr ) ||
268 ( imoread != imo ) ||
272 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chTPFilename );
274 " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
275 iyr, imo , idy ,iyrread, imoread , idyread );
281 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
284 if( chLine[0] !=
'#' && chLine[0] !=
'\n' && chLine[0] !=
'*' )
287 long n = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
294 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
296 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not rewind %s.\n", chTPFilename );
306 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of the transition probability file.\n");
313 lgSentinelReached =
false;
317 if( chLine[0] ==
'*' )
319 lgSentinelReached =
true;
324 if( chLine[0] !=
'#' )
329 if( chLine[0] ==
'\n')
338 ipLo= (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL) - 1;
339 ipHi = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL) - 1;
340 Aul = (double)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
343 fprintf(
ioQQQ,
" PROBLEM End of line reached prematurely in file %s\n",chTPFilename);
344 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
348 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
355 if( (*tr).hasEmis() )
357 fprintf(
ioQQQ,
" PROBLEM duplicate transition found by atmdat_STOUT_readin in %s, "
358 "wavelength=%f\n", chTPFilename,
dBaseStates[intNS][ipHi].energy().WN()
365 (*tr).AddLine2Stack();
366 (*tr).Emis().Aul() = Aul;
367 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
372 fprintf(
ioQQQ,
" PROBLEM File %s contains an invalid line.\n",chTPFilename);
373 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
378 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL );
379 if( !lgSentinelReached )
381 fprintf(
ioQQQ,
" PROBLEM End of data sentinel was not reached in file %s\n",chTPFilename);
382 fprintf(
ioQQQ,
" Check that stars (*****) appear after the last line of data and start in the first column of that line.");
390 printf(
"DEBUG STOUT A's:\n");
394 printf(
"DEBUG.STOUT.A:\t%i\t%i\t%8.2e\n",(*tr).ipLo()+1,(*tr).ipHi()+1,(*tr).Emis().Aul());
404 strcat( chCOLLFilename ,
".coll");
407 ioDATA =
open_data( chCOLLFilename,
"r" );
412 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of %s.\n", chCOLLFilename );
416 iyrread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
417 imoread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
418 idyread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
420 if(( iyrread != iyr ) ||
421 ( imoread != imo ) ||
425 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chCOLLFilename );
427 " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
428 iyr, imo , idy ,iyrread, imoread , idyread );
438 fprintf(
ioQQQ,
" atmdat_STOUT_readin could not read first line of the collision data file.\n");
444 for(
long ipHi=0; ipHi<nMolLevs; ipHi++ )
447 for(
long ipLo=0; ipLo<nMolLevs; ipLo++ )
466 double *temps = NULL;
467 long ipCollider = -1;
468 lgSentinelReached =
false;
474 if( chLine[0] ==
'*' )
476 lgSentinelReached =
true;
481 if( chLine[0] !=
'#' )
486 if( chLine[0] ==
'\n')
490 if(
nMatch(
"TEMP",chLine) )
496 numpoints = (int)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
499 temps = (
double*)
MALLOC(numpoints*
sizeof(
double));
500 memset( temps, 0, (
unsigned long)(numpoints)*
sizeof(
double) );
501 for(
int j = 0; j < numpoints; j++ )
503 temps[j] = (double)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
511 if(
nMatch(
"RATE", chLine) )
514 if(
nMatch(
"ELECTRON",chLine ) )
518 else if(
nMatch(
"PROTON",chLine ) )
524 fprintf(
ioQQQ,
" PROBLEM atmdat_STOUT_readin: Each line of the collision data"
525 "file must specify the collider.\n");
526 fprintf(
ioQQQ,
" Possible colliders are: ELECTRON, PROTON\n");
532 fprintf(
ioQQQ,
" PROBLEM atmdat_STOUT_readin: The collision "
533 "file must specify temperatures before the collision strengths");
538 ipLo= (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL) - 1;
539 ipHi = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL) - 1;
541 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
548 StoutCollData[intNS][ipHi][ipLo][ipCollider].lgIsRate = isRate;
550 StoutCollData[intNS][ipHi][ipLo][ipCollider].ntemps = numpoints;
555 (
double *)
MALLOC((
unsigned long)(numpoints)*
sizeof(
double));
557 (
double *)
MALLOC((
unsigned long)(numpoints)*
sizeof(
double));
559 for(
int j = 0; j < numpoints; j++ )
561 StoutCollData[intNS][ipHi][ipLo][ipCollider].temps[j] = temps[j];
563 StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs[j] = (double)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
569 fprintf(
ioQQQ,
" PROBLEM File %s contains an invalid line.\n",chCOLLFilename);
570 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
576 fprintf(
ioQQQ,
" PROBLEM End of line reached prematurely in file %s\n",chCOLLFilename);
577 fprintf(
ioQQQ,
" The line being read is between the braces {%.*s}\n",
int(strlen(chLine)-1),chLine);
582 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL );
583 if( !lgSentinelReached )
585 fprintf(
ioQQQ,
" PROBLEM End of data sentinel was not reached in file %s\n",chCOLLFilename);
586 fprintf(
ioQQQ,
" Check that stars (*****) appear after the last line of data and start in the first column.");
599 int intsplinepts,intTranType,intxs;
600 long int nMolLevs,nMolExpLevs,nElvlcLines;
601 FILE *ioElecCollData=NULL, *ioProtCollData=NULL;
602 realnum fstatwt,fenergyWN,fWLAng,fenergy,feinsteina;
603 double fScalingParam,fEnergyDiff,*xs,*spl,*spl2;
612 long int *intNewIndex=NULL,*intOldIndex=NULL, *intExperIndex=NULL;
614 bool lgProtonData=
false,lgEneLevOrder;
617 static const int MAX_NUM_LEVELS = 999;
623 strcpy( chEnFilename , chPrefix );
624 strcpy( chTraFilename , chPrefix );
625 strcpy( chEleColFilename , chPrefix );
626 strcpy( chProColFilename , chPrefix );
630 strcat( chEnFilename ,
".elvlc");
633 printf(
"The data file is %s \n",chEnFilename);
637 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chEnFilename);
643 strcat( chTraFilename ,
".wgfa");
646 printf(
"The data file is %s \n",chTraFilename);
650 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chTraFilename);
656 strcat( chEleColFilename ,
".splups");
657 uncaps( chEleColFilename );
659 printf(
"The data file is %s \n",chEleColFilename);
662 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chEleColFilename);
664 ioElecCollData =
open_data( chEleColFilename,
"r" );
667 strcat( chProColFilename ,
".psplups");
668 uncaps( chProColFilename );
670 printf(
"The data file is %s \n",chProColFilename);
673 fprintf(
ioQQQ,
" atmdat_CHIANTI_readin opening %s:",chProColFilename);
676 if( ( ioProtCollData = fopen( chProColFilename ,
"r" ) ) != NULL )
684 lgProtonData =
false;
689 const int eof_col = 5;
691 const int lvl_nrg_col=16;
693 const int lvl_skipto_nrg = 40;
695 const int lvl_eof_to_nrg = lvl_skipto_nrg - eof_col + 1;
698 lgEneLevOrder =
true;
699 if (elvlcstream.is_open())
703 char exptemp[lvl_nrg_col];
704 double tempexpenergy = 0.;
709 elvlcstream.get(otemp,eof_col);
715 elvlcstream.seekg(lvl_eof_to_nrg,ios::cur);
716 elvlcstream.get(exptemp,lvl_nrg_col);
717 tempexpenergy = (
realnum) atof(exptemp);
718 if( tempexpenergy != 0.)
721 elvlcstream.ignore(INT_MAX,
'\n');
724 elvlcstream.seekg(0,ios::beg);
729 printf(
"DEBUG: The file %s contains %li experimental energy levels of the %li total levels. \n",chEnFilename,nMolExpLevs,nElvlcLines);
738 nMolLevs =
MIN3(nMolExpLevs,
atmdat.nChiantiMaxLevelsFe, MAX_NUM_LEVELS );
742 nMolLevs =
MIN3(nMolExpLevs,
atmdat.nChiantiMaxLevels, MAX_NUM_LEVELS );
750 nMolLevs =
MIN3(nElvlcLines,
atmdat.nChiantiMaxLevelsFe,MAX_NUM_LEVELS );
754 nMolLevs =
MIN3(nElvlcLines,
atmdat.nChiantiMaxLevels,MAX_NUM_LEVELS );
760 fprintf(
ioQQQ,
"The number of energy levels is non-positive in datafile %s.\n", chEnFilename );
761 fprintf(
ioQQQ,
"The file must be corrupted.\n" );
762 fclose( ioProtCollData );
769 if(
atmdat.lgChiantiPrint ==
true)
773 fprintf(
ioQQQ,
" Using CHIANTI model %s with %li experimental energy levels of %li available.\n",
778 fprintf(
ioQQQ,
" Using CHIANTI model %s with %li theoretical energy levels of %li available.\n",
791 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
799 for(
long ipHi = 1; ipHi < nMolLevs; ipHi++)
801 for(
long ipLo = 0; ipLo < ipHi; ipLo++)
814 const int lvl_skip_ryd = 15;
824 intExperIndex = (
long int*)
MALLOC((
unsigned long)(nElvlcLines)*
sizeof(
long int));
827 for(
int i = 0; i < nElvlcLines; i++)
829 intExperIndex[i] = -1;
834 for(
long ipLev=0; ipLev<nElvlcLines; ipLev++ )
836 if(elvlcstream.is_open())
838 char thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
839 elvlcstream.seekg(lvl_skipto_nrg,ios::cur);
840 elvlcstream.get(thtemp,lvl_nrg_col);
841 fenergy = (
realnum) atof(thtemp);
852 if( ncounter == nMolLevs)
854 ASSERT( ncounter < nMolLevs );
856 if(fenergy != 0. || ipLev == 0 )
858 dBaseStatesEnergy[ncounter] = fenergy;
859 intExperIndex[ipLev] = ncounter;
864 intExperIndex[ipLev] = -1;
871 if( ipLev == nMolLevs )
874 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
875 elvlcstream.get(obtemp,lvl_nrg_col);
876 if(atof(obtemp) != 0.)
878 fenergy = (
realnum) atof(obtemp);
882 dBaseStatesEnergy[ipLev] = fenergy;
885 elvlcstream.ignore(INT_MAX,
'\n');
890 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEnFilename);
891 fclose( ioProtCollData );
896 for(
long ipLev=1; ipLev<nMolLevs; ipLev++ )
900 if (dBaseStatesEnergy[ipLev] < dBaseStatesEnergy[ipLev-1])
902 lgEneLevOrder =
false;
907 intNewIndex = (
long int*)
MALLOC((
unsigned long)(nElvlcLines)*
sizeof(
long int));
909 if(lgEneLevOrder ==
false)
913 intOldIndex = (
long int*)
MALLOC((
unsigned long)(nMolLevs)*
sizeof(
long int));
915 spsort(dBaseStatesEnergy,nMolLevs,intOldIndex,2,&interror);
917 for(
long i=0; i<nMolLevs; i++ )
919 ASSERT( intOldIndex[i] < nMolLevs );
920 intNewIndex[intOldIndex[i]] = i;
923 if( nMolLevs < nElvlcLines )
930 for(
long i=nMolLevs; i<nElvlcLines; i++)
938 for(
long i=0; i<nMolLevs; i++ )
940 printf(
"The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]);
942 for(
long i=0; i<nMolLevs; i++ )
944 printf(
"The %ld value of the sorted energy array is %f \n",i,dBaseStatesEnergy[i]);
952 for(
long i=0; i<nMolLevs; i++ )
957 if( nMolLevs < nElvlcLines )
961 for(
long i=nMolLevs; i<nElvlcLines; i++)
969 for(
long i=0; i<nMolLevs; i++ )
971 for(
long j=i+1; j<nMolLevs; j++ )
973 ASSERT( intNewIndex[i] != intNewIndex[j] );
983 printf(
"\n\n%s Energy Level Matrix\n",
dBaseSpecies[intNS].chLabel);
984 printf(
"(Original, Experimental, Sorted)\n");
985 for(
long ipLevOld=0; ipLevOld<nElvlcLines; ipLevOld++ )
987 if( intExperIndex[ipLevOld] != -1)
988 printf(
"%li\t%li\t%li\n",ipLevOld+1,intExperIndex[ipLevOld]+1,intNewIndex[intExperIndex[ipLevOld]]+1);
994 elvlcstream.seekg(0,ios::beg);
997 const int lvl_skipto_statwt = 37;
999 const int lvl_statwt_col = 4;
1001 for(
long ipLevOld=0; ipLevOld<nElvlcLines; ipLevOld++ )
1004 if(
atmdat.lgChiantiExp )
1010 if( intExperIndex[ipLevOld] == -1 )
1012 elvlcstream.ignore(INT_MAX,
'\n');
1018 ipLevNew = intNewIndex[intExperIndex[ipLevOld]];
1026 if( intNewIndex[ipLevOld] == -1 )
1028 elvlcstream.ignore(INT_MAX,
'\n');
1033 ipLevNew = intNewIndex[ipLevOld];
1037 char gtemp[lvl_statwt_col],thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
1040 strcpy(
dBaseStates[intNS][ipLevNew].chLabel(),
" ");
1045 elvlcstream.seekg(lvl_skipto_statwt,ios::cur);
1046 elvlcstream.get(gtemp,lvl_statwt_col);
1047 fstatwt = (
realnum)atof(gtemp);
1051 fprintf(
ioQQQ,
" WARNING: A positive non zero value is expected for the "
1052 "statistical weight but was not found in %s"
1053 " level %li\n", chEnFilename,ipLevNew);
1060 elvlcstream.get(thtemp,lvl_nrg_col);
1062 fenergy = (
realnum) atof(thtemp);
1066 if( !
atmdat.lgChiantiExp )
1068 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
1069 elvlcstream.get(obtemp,lvl_nrg_col);
1070 if(atof(obtemp) != 0.)
1072 fenergy = (
realnum) atof(obtemp);
1075 elvlcstream.ignore(INT_MAX,
'\n');
1076 dBaseStates[intNS][ipLevNew].energy().set(fenergy,
"cm^-1");
1079 elvlcstream.close();
1087 int ipHi = (*tr).ipHi();
1088 int ipLo = (*tr).ipLo();
1091 (*tr).EnergyWN() = fenergyWN;
1100 printf(
"\n%s Stored Energy Levels\n",
dBaseSpecies[intNS].chLabel);
1101 printf(
"(Index,Energy in WN, Statistical Weight)\n");
1102 for(
long ipLo = 0; ipLo < nMolLevs; ipLo++ )
1114 if (wgfastream.is_open())
1117 char otemp[eof_col];
1120 wgfastream.get(otemp,eof_col);
1121 wgfastream.ignore(INT_MAX,
'\n');
1125 wgfastream.seekg(0,ios::beg);
1128 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chTraFilename);
1131 const int line_index_col = 6;
1133 const int line_nrg_to_eina =15;
1135 const int line_eina_col = 16;
1136 char lvltemp[line_index_col];
1138 if (wgfastream.is_open())
1140 for (
long i = 0;i<wgfarows;i++)
1142 wgfastream.get(lvltemp,line_index_col);
1143 long ipLo = atoi(lvltemp)-1;
1144 wgfastream.get(lvltemp,line_index_col);
1145 long ipHi = atoi(lvltemp)-1;
1147 if(
atmdat.lgChiantiExp )
1152 if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1154 wgfastream.ignore(INT_MAX,
'\n');
1159 ipLo = intNewIndex[intExperIndex[ipLo]];
1160 ipHi = intNewIndex[intExperIndex[ipHi]];
1168 if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
1170 wgfastream.ignore(INT_MAX,
'\n');
1175 ipLo = intNewIndex[ipLo];
1176 ipHi = intNewIndex[ipHi];
1180 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1183 wgfastream.ignore(INT_MAX,
'\n');
1189 fprintf(
ioQQQ,
" WARNING: Upper level = lower for a radiative transition in %s. Ignoring.\n", chTraFilename );
1190 wgfastream.ignore(INT_MAX,
'\n');
1211 char trantemp[lvl_nrg_col];
1212 wgfastream.get(trantemp,lvl_nrg_col);
1213 fWLAng = (
realnum)atof(trantemp);
1219 if(
atmdat.lgChiantiExp )
1231 long ipLoTemp = ipLo;
1232 long ipHiTemp = ipHi;
1233 ipHi =
max( ipHiTemp, ipLoTemp );
1234 ipLo =
min( ipHiTemp, ipLoTemp );
1235 if(
atmdat.lgChiantiPrint )
1237 fprintf(
ioQQQ,
" WARNING: Swapped level indices for species %6s, indices %3li %3li with wavelength %e \n",
1238 dBaseSpecies[intNS].chLabel, ipLoTemp, ipHiTemp,fWLAng );
1241 else if(
atmdat.lgChiantiPrint )
1243 fprintf(
ioQQQ,
" WARNING: Upper and Lower Indices are reversed.Species: %6s, Indices: %3li %3li Wavelength: %e \n",
1251 wgfastream.ignore(INT_MAX,
'\n');
1252 if(
atmdat.lgChiantiPrint )
1254 fprintf(
ioQQQ,
"WARNING: Skipping Transition %li to %li of %s.\n",ipLo+1,ipHi+1,
dBaseSpecies[intNS].chLabel);
1265 long ipLoTemp = ipLo;
1266 long ipHiTemp = ipHi;
1267 ipHi =
max( ipHiTemp, ipLoTemp );
1268 ipLo =
min( ipHiTemp, ipLoTemp );
1269 fprintf(
ioQQQ,
" WARNING: Swapped level indices for species %6s, indices %3li %3li with wavelength %e \n",
1270 dBaseSpecies[intNS].chLabel, ipLoTemp, ipHiTemp,fWLAng );
1277 if( !
atmdat.lgChiantiExp && fWLAng <= 0. )
1285 wgfastream.seekg(line_nrg_to_eina,ios::cur);
1286 wgfastream.get(trantemp,line_eina_col);
1287 feinsteina = (
realnum)atof(trantemp);
1288 if( feinsteina < 1e-20 )
1290 static bool notPrintedYet =
true;
1291 if( notPrintedYet &&
atmdat.lgChiantiPrint)
1293 fprintf(
ioQQQ,
" WARNING: Radiative rate(s) equal to zero in %s.\n", chTraFilename );
1294 notPrintedYet =
false;
1296 wgfastream.ignore(INT_MAX,
'\n');
1302 wgfastream.getline(chLine,INT_MAX);
1304 if(
nMatch(
"auto", chLine) )
1306 if( (*tr).hasEmis() )
1308 (*tr).Emis().AutoIonizFrac() =
1309 feinsteina/((*tr).Emis().Aul() + feinsteina);
1310 ASSERT( (*tr).Emis().AutoIonizFrac() >= 0. );
1311 ASSERT( (*tr).Emis().AutoIonizFrac() <= 1. );
1316 if( (*tr).hasEmis() )
1318 fprintf(
ioQQQ,
" PROBLEM duplicate transition found by atmdat_chianti in %s, "
1319 "wavelength=%f\n", chTraFilename,fWLAng);
1320 fclose( ioProtCollData );
1323 (*tr).AddLine2Stack();
1324 (*tr).Emis().Aul() = feinsteina;
1325 (*tr).WLAng() = fWLAng;
1327 fenergyWN = (
realnum)(1e+8/fWLAng);
1331 (*tr).EnergyWN() = fenergyWN;
1333 (*tr).Emis().gf() = (
realnum)
GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
1337 fprintf(
ioQQQ,
"The lower level is %ld \n",ipLo);
1338 fprintf(
ioQQQ,
"The upper level is %ld \n",ipHi);
1339 fprintf(
ioQQQ,
"The Einstein A is %f \n",(*tr).Emis().Aul());
1340 fprintf(
ioQQQ,
"The wavelength of the transition is %f \n",(*tr).WLAng() );
1344 else fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chTraFilename);
1349 for(
long ipHi=0; ipHi<nMolLevs; ipHi++ )
1352 for(
long ipLo=0; ipLo<nMolLevs; ipLo++ )
1375 for(
long ipCollider=0; ipCollider < 1; ipCollider++ )
1381 strcpy( chFilename, chEleColFilename );
1387 fprintf(
ioQQQ,
" WARNING: Chianti proton collision data have different format than electron data files!\n" );
1388 strcpy( chFilename, chProColFilename );
1396 fstream splupsstream;
1400 const int cs_eof_col = 4;
1402 const int cs_index_col = 4;
1404 const int cs_trantype_col = 4;
1407 const int cs_values_col = 11;
1409 if (splupsstream.is_open())
1413 long splupslines = -1;
1414 char otemp[cs_eof_col];
1417 splupsstream.get(otemp,cs_eof_col);
1418 splupsstream.ignore(INT_MAX,
'\n');
1422 splupsstream.seekg(0,ios::beg);
1424 for (
int m = 0;m<splupslines;m++)
1428 splupsstream.seekg(6,ios::cur);
1432 splupsstream.get(otemp,cs_index_col);
1433 long ipLo = atoi(otemp)-1;
1434 splupsstream.get(otemp,cs_index_col);
1435 long ipHi = atoi(otemp)-1;
1440 if(
atmdat.lgChiantiExp )
1442 if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1444 splupsstream.ignore(INT_MAX,
'\n');
1449 ipLo = intNewIndex[intExperIndex[ipLo]];
1450 ipHi = intNewIndex[intExperIndex[ipHi]];
1459 if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
1461 splupsstream.ignore(INT_MAX,
'\n');
1466 ipLo = intNewIndex[ipLo];
1467 ipHi = intNewIndex[ipHi];
1471 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1474 splupsstream.ignore(INT_MAX,
'\n');
1478 ASSERT( ipLo < nMolLevs );
1479 ASSERT( ipHi < nMolLevs );
1481 splupsstream.get(otemp,cs_trantype_col);
1482 intTranType = atoi(otemp);
1483 char qtemp[cs_values_col];
1484 splupsstream.get(qtemp,cs_values_col);
1486 splupsstream.get(qtemp,cs_values_col);
1487 fEnergyDiff = atof(qtemp);
1489 splupsstream.get(qtemp,cs_values_col);
1490 fScalingParam = atof(qtemp);
1492 ASSERT( ipLo < nMolLevs );
1493 ASSERT( ipHi < nMolLevs );
1497 const int CHIANTI_SPLINE_MAX=9, CHIANTI_SPLINE_MIN=5;
1502 (
double *)
MALLOC((
unsigned long)(CHIANTI_SPLINE_MAX)*
sizeof(
double));
1504 (
double *)
MALLOC((
unsigned long)(CHIANTI_SPLINE_MAX)*
sizeof(
double));
1507 for( intsplinepts=0; intsplinepts<=CHIANTI_SPLINE_MAX; intsplinepts++ )
1510 char p = splupsstream.peek();
1515 if( (intsplinepts != CHIANTI_SPLINE_MAX) && (intsplinepts != CHIANTI_SPLINE_MIN) )
1517 static bool notPrintedYet =
true;
1520 fprintf(
ioQQQ,
" WARNING: Wrong number of spline points in %s.\n", chFilename);
1521 notPrintedYet =
false;
1523 for(
long j=intsplinepts-1; j<CHIANTI_SPLINE_MAX; j++ )
1531 if( intsplinepts >= CHIANTI_SPLINE_MAX )
1533 fprintf(
ioQQQ,
" WARNING: More spline points than expected in %s, indices %3li %3li. Ignoring extras.\n", chFilename, ipHi, ipLo );
1536 ASSERT( intsplinepts < CHIANTI_SPLINE_MAX );
1539 splupsstream.get(qtemp,cs_values_col);
1545 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intsplinepts] = temp;
1550 ASSERT( (intsplinepts == CHIANTI_SPLINE_MAX) || (intsplinepts == CHIANTI_SPLINE_MIN) );
1559 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].ScalingParam = fScalingParam;
1563 xs = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(
double));
1564 spl = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(
double));
1565 spl2 = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(
double));
1568 if(intsplinepts == CHIANTI_SPLINE_MIN)
1570 for(intxs=0;intxs<CHIANTI_SPLINE_MIN;intxs++)
1572 xs[intxs] = 0.25*intxs;
1573 spl[intxs] =
AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
1576 else if(intsplinepts == CHIANTI_SPLINE_MAX)
1578 for(intxs=0;intxs<CHIANTI_SPLINE_MAX;intxs++)
1580 xs[intxs] = 0.125*intxs;
1581 spl[intxs] =
AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
1587 spline(xs, spl,intsplinepts,2e31,2e31,spl2);
1590 for(
long i=0; i<intsplinepts; i++)
1598 splupsstream.ignore(INT_MAX,
'\n');
1600 splupsstream.close();
1605 free( dBaseStatesEnergy );
1606 free( intNewIndex );
1609 free( intExperIndex );
1613 fclose( ioElecCollData );
1615 fclose( ioProtCollData );