11 #define DEBUGSTATE false
17 int intCS,intCollIndex = -10000,intsplinepts,intTranType,intxs;
19 long int i,j,nMolLevs,intCollTran,intcollindex,ipLo,ipHi;
20 FILE *atmolLevDATA , *atmolTraDATA=NULL,*atmolEleColDATA=NULL,*atmolProColDATA=NULL;
21 realnum fstatwt,fenergyK,fenergyWN,fWLAng,fenergy,feinsteina;
22 double fScalingParam,fEnergyDiff,fGF,*xs,*spl,*spl2;
30 long ECS,SWS,LLS,ULS,EAS,WOTS;
40 long int *intNewIndex=NULL,*intOldIndex=NULL;
42 bool lgProtonData=
false,lgEneLevOrder;
55 strcpy( chEnFilename ,
Species[intNS].chptrSpName );
56 strcat( chEnFilename ,
".elvlc");
59 printf(
"The data file is %s \n",chEnFilename);
63 fprintf(
ioQQQ,
" moldat_readin opening %s:",chEnFilename);
65 atmolLevDATA =
open_data( chEnFilename,
"r" );
68 strcpy( chTraFilename ,
Species[intNS].chptrSpName );
69 strcat( chTraFilename ,
".wgfa");
72 printf(
"The data file is %s \n",chTraFilename);
76 fprintf(
ioQQQ,
" moldat_readin opening %s:",chTraFilename);
78 atmolTraDATA =
open_data( chTraFilename,
"r" );
81 strcpy( chEleColFilename ,
Species[intNS].chptrSpName );
82 strcat( chEleColFilename ,
".splups");
83 uncaps( chEleColFilename );
85 printf(
"The data file is %s \n",chEleColFilename);
88 fprintf(
ioQQQ,
" moldat_readin opening %s:",chEleColFilename);
90 atmolEleColDATA =
open_data( chEleColFilename,
"r" );
93 strcpy( chProColFilename ,
Species[intNS].chptrSpName );
94 strcat( chProColFilename ,
".psplups");
95 uncaps( chProColFilename );
97 printf(
"The data file is %s \n",chProColFilename);
100 fprintf(
ioQQQ,
" moldat_readin opening %s:",chProColFilename);
102 if( ( atmolProColDATA = fopen( chProColFilename ,
"r" ) ) != NULL )
105 fclose( atmolProColDATA );
106 atmolProColDATA = NULL;
110 lgProtonData =
false;
114 lgEneLevOrder =
true;
115 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , atmolLevDATA ) != NULL )
117 chDesired = strtok(chLine,
" \t");
118 intCS = atoi(chDesired);
128 printf(
"The number of energy levels is %li",nMolLevs);
132 fprintf(
ioQQQ,
"The number of energy levels is non-positive in datafile %s.\n", chEnFilename );
133 fprintf(
ioQQQ,
"The file must be corrupted.\n" );
146 for( ipHi = 1; ipHi < nMolLevs; ipHi++)
149 for(ipLo = 0; ipLo < ipHi; ipLo++)
157 for( intcollindex = 0; intcollindex<
NUM_COLLIDERS; intcollindex++ )
163 for( ipHi = 0; ipHi<nMolLevs; ipHi++ )
166 for( ipLo = 0; ipLo<nMolLevs; ipLo++)
175 for( ipHi = 0; ipHi<nMolLevs; ipHi++ )
178 for( ipLo = 0; ipLo<nMolLevs; ipLo++)
185 if( fseek( atmolLevDATA , 0 , SEEK_SET ) != 0 )
187 fprintf(
ioQQQ,
" moldat_readin could not rewind %s.\n",chEnFilename);
193 for( i=0; i<nMolLevs; i++ )
195 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , atmolLevDATA ) == NULL )
197 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEnFilename);
200 fenergy = (
realnum) atof(&chLine[ECS]);
201 atmolStatesEnergy[i] = fenergy;
207 if (atmolStatesEnergy[i] < atmolStatesEnergy[i-1])
209 lgEneLevOrder =
false;
216 for(i=0;i<nMolLevs;i++)
218 printf(
"The %ld value of the unsorted energy array is %f \n",i,atmolStatesEnergy[i]);
222 intNewIndex = (
long int*)
MALLOC((
unsigned long)(nMolLevs)*
sizeof(
long int));
224 if(lgEneLevOrder ==
false)
228 intOldIndex = (
long int*)
MALLOC((
unsigned long)(nMolLevs)*
sizeof(
long int));
230 spsort(atmolStatesEnergy,nMolLevs,intOldIndex,2,&interror);
232 for( i=0; i<nMolLevs; i++ )
234 ASSERT( intOldIndex[i] < nMolLevs );
235 intNewIndex[intOldIndex[i]] = i;
239 for(i=0;i<nMolLevs;i++)
241 printf(
"The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]);
243 for( i=0; i<nMolLevs; i++ )
245 printf(
"The %ld value of the sorted energy array is %f \n",i,atmolStatesEnergy[i]);
254 for( i=0; i<nMolLevs; i++ )
261 for( i=0; i<nMolLevs; i++ )
263 for( j=i+1; j<nMolLevs; j++ )
265 ASSERT( intNewIndex[i] != intNewIndex[j] );
271 if( fseek( atmolLevDATA , 0 , SEEK_SET ) != 0 )
273 fprintf(
ioQQQ,
" moldat_readin could not rewind %s.\n",chEnFilename);
277 for( i=0; i<nMolLevs; i++ )
279 long j = intNewIndex[i];
281 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , atmolLevDATA ) == NULL )
283 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEnFilename);
291 fstatwt = (
realnum)atof(&chLine[SWS]);
294 fprintf(
ioQQQ,
" A positive non zero value is expected for the statistical weight .\n");
295 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEnFilename);
300 fenergy = (
realnum) atof(&chLine[ECS]);
305 fprintf(
ioQQQ,
"The %ld value of g after populating is %f \n",
307 fprintf(
ioQQQ,
"The %ld value of energy after populating is %f \n",
314 for( i=1; i<nMolLevs; i++ )
331 if(
read_whole_line( chLine , (
int)
sizeof(chLine) ,atmolTraDATA ) == NULL )
333 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chTraFilename);
337 while(atoi(chLine)!= -1)
339 ipLo = atoi(&chLine[LLS])-1;
340 ipHi = atoi(&chLine[ULS])-1;
343 ipLo = intNewIndex[ipLo];
344 ipHi = intNewIndex[ipHi];
348 fWLAng = (
realnum)atof(&chLine[WOTS]);
359 feinsteina = (
realnum)atof(&chLine[EAS]);
363 fenergyWN = (
realnum)(1e+8/fWLAng);
375 fprintf(
ioQQQ,
"The lower level is %ld \n",ipLo);
376 fprintf(
ioQQQ,
"The upper level is %ld \n",ipHi);
377 fprintf(
ioQQQ,
"The Einstein A is %f \n",
atmolTrans[intNS][ipHi][ipLo].Emis->Aul);
378 fprintf(
ioQQQ,
"The wavelength of the transition is %f \n",
atmolTrans[intNS][ipHi][ipLo].WLAng );
381 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , atmolTraDATA ) == NULL )
383 fprintf(
ioQQQ,
" Data is expected on this line.\n");
384 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chTraFilename);
392 for( i=0; i<nMolLevs; i++ )
395 for( j=0; j<nMolLevs; j++ )
421 if(
read_whole_line( chLine , (
int)
sizeof(chLine) ,atmolEleColDATA ) == NULL )
423 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEleColFilename);
427 while(atoi(chLine)!= -1)
430 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , atmolEleColDATA ) == NULL )
432 fprintf(
ioQQQ,
" Data is expected on this line.\n");
433 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEleColFilename);
440 printf(
"\n The number of collisional transitions is %li",intCollTran);
443 if( fseek( atmolEleColDATA , 0 , SEEK_SET ) != 0 )
445 fprintf(
ioQQQ,
" moldat_readin could not rewind %s.\n",chEleColFilename);
451 if(
read_whole_line( chLine , (
int)
sizeof(chLine) ,atmolEleColDATA ) == NULL )
453 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEleColFilename);
458 if(atoi(chLine)== -1)
460 fprintf(
ioQQQ,
" If the file %s exists,the first line cannot be -1 .\n",chEleColFilename);
464 while(atoi(chLine)!= -1)
469 (void)
FFmtRead( chLine, &i,
sizeof(chLine), &lgEOL );
470 (void)
FFmtRead( chLine, &i,
sizeof(chLine), &lgEOL );
472 ipLo = (long)
FFmtRead( chLine, &i,
sizeof(chLine), &lgEOL ) - 1;
473 ipHi = (long)
FFmtRead( chLine, &i,
sizeof(chLine), &lgEOL ) - 1;
474 intTranType = (long)
FFmtRead( chLine, &i,
sizeof(chLine), &lgEOL );
480 fScalingParam = (
realnum)
FFmtRead( chLine, &i,
sizeof(chLine), &lgEOL );
483 ipLo = intNewIndex[ipLo];
484 ipHi = intNewIndex[ipHi];
486 ASSERT( ipLo < nMolLevs );
487 ASSERT( ipHi < nMolLevs );
494 (
double *)
MALLOC((
unsigned long)(9)*
sizeof(
double));
496 (
double *)
MALLOC((
unsigned long)(9)*
sizeof(
double));
499 for( intsplinepts=0; intsplinepts<10; intsplinepts++ )
501 double temp = (double)
FFmtRead( chLine, &i,
sizeof(chLine), &lgEOL );
507 if( (intsplinepts != 5) && (intsplinepts != 9) )
509 fprintf(
ioQQQ,
" There should have been 5 or 9 spline points. Sorry.\n" );
510 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
518 ASSERT( intsplinepts < 9 );
523 ASSERT( (intsplinepts == 5) || (intsplinepts == 9) );
538 xs = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
539 spl = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
540 spl2 = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
542 if(intsplinepts == 5)
544 for(intxs=0;intxs<5;intxs++)
546 xs[intxs] = 0.25*intxs;
550 else if(intsplinepts == 9)
552 for(intxs=0;intxs<9;intxs++)
554 xs[intxs] = 0.125*intxs;
561 spline(xs, spl,intsplinepts,2e31,2e31,spl2);
564 for(i=0;i<intsplinepts;i++)
573 if(
read_whole_line( chLine , (
int)
sizeof(chLine) ,atmolEleColDATA ) == NULL )
575 fprintf(
ioQQQ,
" The data file %s is corrupted .\n",chEleColFilename);
591 fprintf(
ioQQQ,
"Collisional Data: %li\t%li\t%f\t%f\t%f",
601 fprintf(
ioQQQ,
"\n" );
604 free( atmolStatesEnergy );
607 fclose( atmolLevDATA );
608 fclose( atmolTraDATA );
609 fclose( atmolEleColDATA );