cloudy trunk
Loading...
Searching...
No Matches
atmdat_chianti.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3#include "cddefines.h"
4#include "lines_service.h"
5#include "physconst.h"
6#include "taulines.h"
7#include "trace.h"
8#include "string.h"
9#include "thirdparty.h"
10#include "rfield.h"
11#include "atmdat.h"
12
13static const bool DEBUGSTATE = false;
14// minimum energy difference (wavenumbers) we will accept
15const double ENERGY_MIN_WN = 1e-10;
16
17void atmdat_STOUT_readin( long intNS, char *chPrefix )
18{
19 DEBUG_ENTRY( "atmdat_STOUT_readin()" );
20
21 long int nMolLevs;
22 char chLine[FILENAME_PATH_LENGTH_2],
23 chNRGFilename[FILENAME_PATH_LENGTH_2],
24 chTPFilename[FILENAME_PATH_LENGTH_2],
25 chCOLLFilename[FILENAME_PATH_LENGTH_2];
26
27 static const int MAX_NUM_LEVELS = 999;
28
29 dBaseSpecies[intNS].lgMolecular = false;
30 dBaseSpecies[intNS].lgLTE = false;
31 dBaseSpecies[intNS].lgLAMDA = false;
32
33 strcpy( chNRGFilename , chPrefix );
34 strcpy( chTPFilename , chPrefix );
35 strcpy( chCOLLFilename , chPrefix );
36
37 /*Open the energy levels file*/
38 strcat( chNRGFilename , ".nrg");
39 uncaps( chNRGFilename );
40 if(DEBUGSTATE)
41 printf("The data file is %s \n",chNRGFilename);
42
43 /*Open the files*/
44 if( trace.lgTrace )
45 fprintf( ioQQQ," atmdat_STOUT_readin opening %s:",chNRGFilename);
46
47 FILE *ioDATA;
48 ioDATA = open_data( chNRGFilename, "r" );
49 long int i = 0;
50 bool lgEOL = false;
51 long index = 0;
52 double nrg = 0.0;
53 double stwt = 0.0;
54
55 /* first line is a version number - now confirm that it is valid */
56 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
57 {
58 fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chNRGFilename );
60 }
61 i = 1;
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);
67
68 if(( iyrread != iyr ) ||
69 ( imoread != imo ) ||
70 ( idyread != idy ) )
71 {
72 fprintf( ioQQQ,
73 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chNRGFilename );
74 fprintf( ioQQQ,
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 );
78 }
79
80 nMolLevs = 0;
81 //Count number of levels
82 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
83 {
84 /* # - comment, *** ends data */
85 if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
86 {
87 i = 1;
88 long n = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
89 if( n < 0 )
90 break;
91 nMolLevs++;
92 }
93 }
94
95 /* now rewind the file so we can read it a second time*/
96 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
97 {
98 fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chNRGFilename );
100 }
101 //Skip the magic numbers this time
102 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
103
104 long HighestIndexInFile = nMolLevs;
105
106 nMolLevs = MIN3(atmdat.nStoutMaxLevels,nMolLevs, MAX_NUM_LEVELS );
107
108 if( atmdat.lgStoutPrint == true)
109 {
110 fprintf( ioQQQ," Using STOUT model %s with %li levels of %li available.\n",
111 dBaseSpecies[intNS].chLabel , nMolLevs , HighestIndexInFile );
112 }
113
114 dBaseSpecies[intNS].numLevels_max = nMolLevs;
115 dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
116
117 /*Resize the States array*/
118 dBaseStates[intNS].resize(nMolLevs);
119 /*Zero out the maximum wavenumber for each species */
120 dBaseSpecies[intNS].maxWN = 0.;
121
122 /* allocate the Transition array*/
123 ipdBaseTrans[intNS].reserve(nMolLevs);
124 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
125 ipdBaseTrans[intNS].reserve(ipHi,ipHi);
126 ipdBaseTrans[intNS].alloc();
127 dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
128 dBaseTrans[intNS].states() = &dBaseStates[intNS];
129 dBaseTrans[intNS].chLabel() = "Stout";
130
131 //This is creating transitions that we don't have collision data for
132 //Maybe use gbar or keep all of the Fe 2 even if it was assumed (1e-5)
133 int ndBase = 0;
134 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
135 {
136 for( long ipLo = 0; ipLo < ipHi; ipLo++)
137 {
138 ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
139 dBaseTrans[intNS][ndBase].Junk();
140 dBaseTrans[intNS][ndBase].setLo(ipLo);
141 dBaseTrans[intNS][ndBase].setHi(ipHi);
142 ++ndBase;
143 }
144 }
145
146
147 /* Check for an end of file sentinel */
148 bool lgSentinelReached = false;
149
150 //Read the first line of data
151 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
152 {
153 fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the energy level file.\n");
155 }
156 //Read the remaining lines of the energy level file
157 do
158 {
159 i = 1;
160
161 /* Stop on *** */
162 if( chLine[0] == '*' )
163 {
164 lgSentinelReached = true;
165 break;
166 }
167 //Comments start with #, skip them
168 if( chLine[0] != '#' )
169 {
170 /* Skip blank lines */
171 if( chLine[0] == '\n')
172 continue;
173
174 index = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) -1 ;
175
176 if( index < 0 )
177 {
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);
181 }
182
183 nrg = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
184 stwt = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
185
186 if( lgEOL )
187 {
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);
191 }
192
193 if( index >= nMolLevs )
194 {
195 //Skip unused levels
196 continue;
197 }
198
199 dBaseStates[intNS][index].energy().set(nrg,"cm^-1");
200 dBaseStates[intNS][index].g() = stwt;
201
202 /*information needed for label*/
203 strcpy(dBaseStates[intNS][index].chLabel(), " ");
204 strncpy(dBaseStates[intNS][index].chLabel(),dBaseSpecies[intNS].chLabel, 4);
205 dBaseStates[intNS][index].chLabel()[4] = '\0';
206 }
207 }
208 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
209 if( !lgSentinelReached )
210 {
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");
214 }
215 fclose(ioDATA);
216
217 if( DEBUGSTATE )
218 {
219 /*Test whether energy levels are read in properly*/
220 printf("DEBUG STOUT ENERGY LEVELS:\n");
221 for( int i = 0; i< nMolLevs; i++ )
222 {
223 printf("DEBUG\t%i\t%11.4e\t%3.1f\n",i+1,dBaseStates[intNS][i].energy().WN(),dBaseStates[intNS][i].g());
224 }
225 }
226
227 /* fill in all transition energies, can later overwrite for specific radiative transitions */
228 double fenergyWN = 0.;
229 for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
230 tr!= dBaseTrans[intNS].end(); ++tr)
231 {
232 int ipHi = (*tr).ipHi();
233 int ipLo = (*tr).ipLo();
234 ASSERT(ipHi > ipLo );
235 fenergyWN = dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN();
236 if( fenergyWN <= 0.)
237 {
238 printf("\nWARNING: The %s transition between level %i and %i has zero or negative energy.\n",
239 dBaseStates[intNS][ipHi].chLabel(),ipLo+1,ipHi+1);
240 printf("Check the Stout energy level data file (*.nrg) of this species for missing or duplicate energies.\n");
241 //cdEXIT(EXIT_FAILURE);
242 }
243 (*tr).EnergyWN() = (realnum)MAX2(ENERGY_MIN_WN,fenergyWN);
244 (*tr).WLAng() = (realnum)(1e+8/(*tr).EnergyWN()/RefIndex((*tr).EnergyWN()));
245 dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,(*tr).EnergyWN());
246 }
247
248 /******************************************************
249 ************* Transition Probability File ************
250 ******************************************************/
251 strcat( chTPFilename , ".tp");
252 uncaps( chTPFilename );
253
254 ioDATA = open_data( chTPFilename, "r" );
255
256 /* first line is a version number - now confirm that it is valid */
257 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
258 {
259 fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chTPFilename );
261 }
262 i = 1;
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);
266
267 if(( iyrread != iyr ) ||
268 ( imoread != imo ) ||
269 ( idyread != idy ) )
270 {
271 fprintf( ioQQQ,
272 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chTPFilename );
273 fprintf( ioQQQ,
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 );
277 }
278
279 long numtrans = 0;
280 //Count number of transitions
281 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
282 {
283 /* # is comment, *** is end of data */
284 if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
285 {
286 i = 1;
287 long n = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
288 if( n < 0 )
289 break;
290 numtrans++;
291 }
292 }
293 /* now rewind the file so we can read it a second time*/
294 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
295 {
296 fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chTPFilename );
298 }
299 //Skip the magic numbers this time
300 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
301
302
303 //Read the first line of data
304 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
305 {
306 fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the transition probability file.\n");
308 }
309
310 long ipLo = 0;
311 long ipHi = 0;
312 double Aul = 0.0;
313 lgSentinelReached = false;
314 //Read the remaining lines of the transition probability file
315 do
316 {
317 if( chLine[0] == '*' )
318 {
319 lgSentinelReached = true;
320 break;
321 }
322
323 //Comments start with #, skip them
324 if( chLine[0] != '#' )
325 {
326 i = 1;
327
328 /* skip null lines */
329 if( chLine[0] == '\n')
330 continue;
331
332 //This means last data column has Aul.
333 if( nMatch("A",chLine) )
334 {
335 /* reset read pointer */
336 i = 1;
337
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);
341 if( lgEOL )
342 {
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);
346 }
347
348 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
349 {
350 // skip these lines
351 continue;
352 }
353
354 TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
355 if( (*tr).hasEmis() )
356 {
357 fprintf(ioQQQ," PROBLEM duplicate transition found by atmdat_STOUT_readin in %s, "
358 "wavelength=%f\n", chTPFilename,dBaseStates[intNS][ipHi].energy().WN()
359 - dBaseStates[intNS][ipLo].energy().WN());
361 }
362
363 if( (*tr).EnergyWN() > ENERGY_MIN_WN )
364 {
365 (*tr).AddLine2Stack();
366 (*tr).Emis().Aul() = Aul;
367 (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
368 }
369 }
370 else
371 {
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);
375 }
376 }
377 }
378 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
379 if( !lgSentinelReached )
380 {
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.");
384 }
385 fclose(ioDATA);
386
387 if( DEBUGSTATE )
388 {
389 /*Test whether stout A's are read in properly */
390 printf("DEBUG STOUT A's:\n");
391 for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
392 tr!= dBaseTrans[intNS].end(); ++tr)
393 {
394 printf("DEBUG.STOUT.A:\t%i\t%i\t%8.2e\n",(*tr).ipLo()+1,(*tr).ipHi()+1,(*tr).Emis().Aul());
395 }
396 }
397
398
399
400 /******************************************************
401 ************* Collision Data File ********************
402 ******************************************************/
403
404 strcat( chCOLLFilename , ".coll");
405 uncaps( chCOLLFilename );
406
407 ioDATA = open_data( chCOLLFilename, "r" );
408
409 /* first line is a version number - now confirm that it is valid */
410 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
411 {
412 fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chCOLLFilename );
414 }
415 i = 1;
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);
419
420 if(( iyrread != iyr ) ||
421 ( imoread != imo ) ||
422 ( idyread != idy ) )
423 {
424 fprintf( ioQQQ,
425 " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chCOLLFilename );
426 fprintf( ioQQQ,
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 );
430 }
431
432 /****** Could add ability to count number of temperature changes, electron CS, and proton CS ****/
433
434
435 //Read the first line of data
436 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
437 {
438 fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the collision data file.\n");
440 }
441
442 /* Malloc space for collision strengths */
443 StoutCollData[intNS] = (StoutColls***)MALLOC(nMolLevs *sizeof(StoutColls**));
444 for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
445 {
446 StoutCollData[intNS][ipHi] = (StoutColls **)MALLOC((unsigned long)(nMolLevs)*sizeof(StoutColls *));
447 for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
448 {
449 StoutCollData[intNS][ipHi][ipLo] =
450 (StoutColls *)MALLOC((unsigned long)(ipNCOLLIDER)*sizeof(StoutColls ));
451
452 for( long k=0; k<ipNCOLLIDER; k++ )
453 {
454 /* initialize all spline variables */
455 StoutCollData[intNS][ipHi][ipLo][k].ntemps = -1;
456 StoutCollData[intNS][ipHi][ipLo][k].temps = NULL;
457 StoutCollData[intNS][ipHi][ipLo][k].collstrs = NULL;
458 StoutCollData[intNS][ipHi][ipLo][k].lgIsRate = false;
459 }
460 }
461 }
462
463 ipLo = 0;
464 ipHi = 0;
465 int numpoints = 0;
466 double *temps = NULL;
467 long ipCollider = -1;
468 lgSentinelReached = false;
469
470 //Read the remaining lines of the collision data file
471 do
472 {
473 /* Stop on *** */
474 if( chLine[0] == '*' )
475 {
476 lgSentinelReached = true;
477 break;
478 }
479
480 //Comments start with #, skip them
481 if( chLine[0] != '#' )
482 {
483 i = 1;
484
485 /* Skip blank lines */
486 if( chLine[0] == '\n')
487 continue;
488
489 //This is a temperature line
490 if( nMatch("TEMP",chLine) )
491 {
492 if( temps != NULL)
493 free( temps );
494
495 i = 1;
496 numpoints = (int)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
497 ASSERT( numpoints > 0 );
498
499 temps = (double*)MALLOC(numpoints*sizeof(double));
500 memset( temps, 0, (unsigned long)(numpoints)*sizeof(double) );
501 for( int j = 0; j < numpoints; j++ )
502 {
503 temps[j] = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
504 ASSERT( temps[j] > 0 );
505 }
506 }
507 else if( nMatch("CS",chLine) || nMatch("RATE",chLine) )
508 {
509
510 bool isRate = false;
511 if( nMatch("RATE", chLine) )
512 isRate = true;
513
514 if( nMatch( "ELECTRON",chLine ) )
515 {
516 ipCollider = ipELECTRON;
517 }
518 else if( nMatch( "PROTON",chLine ) )
519 {
520 ipCollider = ipPROTON;
521 }
522 else
523 {
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");
528 }
529
530 if( temps == NULL )
531 {
532 fprintf( ioQQQ, " PROBLEM atmdat_STOUT_readin: The collision "
533 "file must specify temperatures before the collision strengths");
535 }
536
537 i = 1;
538 ipLo= (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
539 ipHi = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
540
541 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
542 {
543 // skip these lines
544 continue;
545 }
546
547 /* Set this as a collision strength not a collision rate */
548 StoutCollData[intNS][ipHi][ipLo][ipCollider].lgIsRate = isRate;
549
550 StoutCollData[intNS][ipHi][ipLo][ipCollider].ntemps = numpoints;
551 ASSERT( numpoints > 0 );
552
553 /*malloc the space here for the temps and collision strengths*/
554 StoutCollData[intNS][ipHi][ipLo][ipCollider].temps =
555 (double *)MALLOC((unsigned long)(numpoints)*sizeof(double));
556 StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs =
557 (double *)MALLOC((unsigned long)(numpoints)*sizeof(double));
558 /* Loop over all but one CS value */
559 for( int j = 0; j < numpoints; j++ )
560 {
561 StoutCollData[intNS][ipHi][ipLo][ipCollider].temps[j] = temps[j];
562 ASSERT( StoutCollData[intNS][ipHi][ipLo][ipCollider].temps[j] > 0 );
563 StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs[j] = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
564 ASSERT( StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs[j] > 0 );
565 }
566 }
567 else
568 {
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);
572 }
573
574 if( lgEOL )
575 {
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);
579 }
580 }
581 }
582 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
583 if( !lgSentinelReached )
584 {
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.");
588 }
589 free(temps);
590 fclose(ioDATA);
591
592 return;
593}
594
595void atmdat_CHIANTI_readin( long intNS, char *chPrefix )
596{
597 DEBUG_ENTRY( "atmdat_CHIANTI_readin()" );
598
599 int intsplinepts,intTranType,intxs;
600 long int nMolLevs,nMolExpLevs,nElvlcLines;// number of experimental and total levels
601 FILE *ioElecCollData=NULL, *ioProtCollData=NULL;
602 realnum fstatwt,fenergyWN,fWLAng,fenergy,feinsteina;
603 double fScalingParam,fEnergyDiff,*xs,*spl,*spl2;
604
605 char chLine[FILENAME_PATH_LENGTH_2] ,
606 chEnFilename[FILENAME_PATH_LENGTH_2],
607 chTraFilename[FILENAME_PATH_LENGTH_2],
608 chEleColFilename[FILENAME_PATH_LENGTH_2],
609 chProColFilename[FILENAME_PATH_LENGTH_2];
610
611 realnum *dBaseStatesEnergy;
612 long int *intNewIndex=NULL,*intOldIndex=NULL, *intExperIndex=NULL;
613 int interror;
614 bool lgProtonData=false,lgEneLevOrder;
615
616 // this is the largest number of levels allowed by the chianti format, I3
617 static const int MAX_NUM_LEVELS = 999;
618
619 dBaseSpecies[intNS].lgMolecular = false;
620 dBaseSpecies[intNS].lgLAMDA = false;
621 dBaseSpecies[intNS].lgLTE = false;
622
623 strcpy( chEnFilename , chPrefix );
624 strcpy( chTraFilename , chPrefix );
625 strcpy( chEleColFilename , chPrefix );
626 strcpy( chProColFilename , chPrefix );
627
628 /*For the CHIANTI DATABASE*/
629 /*Open the energy levels file*/
630 strcat( chEnFilename , ".elvlc");
631 uncaps( chEnFilename );
632 if(DEBUGSTATE)
633 printf("The data file is %s \n",chEnFilename);
634
635 /*Open the files*/
636 if( trace.lgTrace )
637 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEnFilename);
638
639 fstream elvlcstream;
640 open_data(elvlcstream, chEnFilename,mode_r);
641
642 /*Open the transition probabilities file*/
643 strcat( chTraFilename , ".wgfa");
644 uncaps( chTraFilename );
645 if(DEBUGSTATE)
646 printf("The data file is %s \n",chTraFilename);
647
648 /*Open the files*/
649 if( trace.lgTrace )
650 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chTraFilename);
651
652 fstream wgfastream;
653 open_data(wgfastream, chTraFilename,mode_r);
654
655 /*Open the electron collision data*/
656 strcat( chEleColFilename , ".splups");
657 uncaps( chEleColFilename );
658 if(DEBUGSTATE)
659 printf("The data file is %s \n",chEleColFilename);
660 /*Open the files*/
661 if( trace.lgTrace )
662 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEleColFilename);
663
664 ioElecCollData = open_data( chEleColFilename, "r" );
665
666 /*Open the proton collision data*/
667 strcat( chProColFilename , ".psplups");
668 uncaps( chProColFilename );
669 if(DEBUGSTATE)
670 printf("The data file is %s \n",chProColFilename);
671 /*Open the files*/
672 if( trace.lgTrace )
673 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chProColFilename);
674
675 /*We will set a flag here to indicate if the proton collision strengths are available */
676 if( ( ioProtCollData = fopen( chProColFilename , "r" ) ) != NULL )
677 {
678 lgProtonData = true;
679 //fclose( ioProtCollData );
680 //ioProtCollData = NULL;
681 }
682 else
683 {
684 lgProtonData = false;
685 }
686
687 /*Loop finds how many theoretical and experimental levels are in the elvlc file */
688 //eof_col is used get the first 4 charcters per line to find end of file
689 const int eof_col = 5;
690 //length (+1) of the nrg in the elvlc file
691 const int lvl_nrg_col=16;
692 //# of columns skipped from the left to get to nrg start
693 const int lvl_skipto_nrg = 40;
694 /* # of columns to skip from eof check to nrg start */
695 const int lvl_eof_to_nrg = lvl_skipto_nrg - eof_col + 1;
696 nElvlcLines = 0;
697 nMolExpLevs = 1;
698 lgEneLevOrder = true;
699 if (elvlcstream.is_open())
700 {
701 int nj = 0;
702 char otemp[eof_col];
703 char exptemp[lvl_nrg_col];
704 double tempexpenergy = 0.;
705 /*This loop counts the number of valid rows within the elvlc file
706 as well as the number of experimental energy levels.*/
707 while(nj != -1)
708 {
709 elvlcstream.get(otemp,eof_col);
710 nj = atoi(otemp);
711 if( nj == -1)
712 break;
713 nElvlcLines++;
714
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.)
719 nMolExpLevs++;
720
721 elvlcstream.ignore(INT_MAX,'\n');
722
723 }
724 elvlcstream.seekg(0,ios::beg);
725 }
726
727 if(DEBUGSTATE)
728 {
729 printf("DEBUG: The file %s contains %li experimental energy levels of the %li total levels. \n",chEnFilename,nMolExpLevs,nElvlcLines);
730 }
731
732 /* The total number of levels depends on the experimental Chianti switch */
733 if( atmdat.lgChiantiExp )
734 {
735 if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
736 {
737 // Fe is special case with more levels
738 nMolLevs = MIN3(nMolExpLevs, atmdat.nChiantiMaxLevelsFe, MAX_NUM_LEVELS );
739 }
740 else
741 {
742 nMolLevs = MIN3(nMolExpLevs, atmdat.nChiantiMaxLevels, MAX_NUM_LEVELS );
743 }
744 }
745 else
746 {
747 if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
748 {
749 // Fe is special case with more levels
750 nMolLevs = MIN3(nElvlcLines, atmdat.nChiantiMaxLevelsFe,MAX_NUM_LEVELS );
751 }
752 else
753 {
754 nMolLevs = MIN3(nElvlcLines, atmdat.nChiantiMaxLevels,MAX_NUM_LEVELS );
755 }
756 }
757
758 if( nMolLevs <= 0 )
759 {
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 );
764 }
765
766 dBaseSpecies[intNS].numLevels_max = nMolLevs;
767 dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
768
769 if( atmdat.lgChiantiPrint == true)
770 {
771 if( atmdat.lgChiantiExp )
772 {
773 fprintf( ioQQQ," Using CHIANTI model %s with %li experimental energy levels of %li available.\n",
774 dBaseSpecies[intNS].chLabel , nMolLevs , nMolExpLevs );
775 }
776 else
777 {
778 fprintf( ioQQQ," Using CHIANTI model %s with %li theoretical energy levels of %li available.\n",
779 dBaseSpecies[intNS].chLabel , nMolLevs , nElvlcLines );
780 }
781 }
782
783 /*Malloc the energy array to check if the energy levels are in order*/
784 dBaseStatesEnergy = (realnum*)MALLOC((unsigned long)(nMolLevs)*sizeof(realnum));
785
786 /*malloc the States array*/
787 dBaseStates[intNS].resize(nMolLevs);
788
789 /* allocate the Transition array*/
790 ipdBaseTrans[intNS].reserve(nMolLevs);
791 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
792 ipdBaseTrans[intNS].reserve(ipHi,ipHi);
793 ipdBaseTrans[intNS].alloc();
794 dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
795 dBaseTrans[intNS].states() = &dBaseStates[intNS];
796 dBaseTrans[intNS].chLabel() = "Chianti";
797
798 int ndBase = 0;
799 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
800 {
801 for( long ipLo = 0; ipLo < ipHi; ipLo++)
802 {
803 ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
804 dBaseTrans[intNS][ndBase].Junk();
805 dBaseTrans[intNS][ndBase].setLo(ipLo);
806 dBaseTrans[intNS][ndBase].setHi(ipHi);
807 ++ndBase;
808 }
809 }
810
811 /*Reading in the energy and checking that they are in order*/
812 // C++ io works in terms of cursor movement rather than position on line
813 //# of columns to skip over the rydberg energy, we don't use it
814 const int lvl_skip_ryd = 15;
815
816 /*Keep track of which levels have experimental data and then create a vector
817 which relates their indices to the default chianti energy indices.
818 */
819 long ncounter = 0;
820
821 if( atmdat.lgChiantiExp )
822 {
823 //Relate Chianti level indices to a set that only include experimental levels
824 intExperIndex = (long int*)MALLOC((unsigned long)(nElvlcLines)* sizeof(long int));
825
826 //Fill relational vector with bad values
827 for(int i = 0; i < nElvlcLines; i++)
828 {
829 intExperIndex[i] = -1;
830 }
831 }
832
833 //Read in nrg levels to see if they are in order
834 for( long ipLev=0; ipLev<nElvlcLines; ipLev++ )
835 {
836 if(elvlcstream.is_open())
837 {
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);
842
843 if( atmdat.lgChiantiExp )
844 {
845 /* Go through the entire level list selectively choosing only experimental level energies.
846 * Store them, not zeroes, in order using ncounter to count the index.
847 * Any row on the level list where there is no experimental energy, put a -1 in the relational vector.
848 * If it is a valid experimental energy level store the new ncounter index.
849 */
850
851 //Once we collect enough experimental levels stop looking for more.
852 if( ncounter == nMolLevs)
853 break;
854 ASSERT( ncounter < nMolLevs );
855
856 if(fenergy != 0. || ipLev == 0 )
857 {
858 dBaseStatesEnergy[ncounter] = fenergy;
859 intExperIndex[ipLev] = ncounter;
860 ncounter++;
861 }
862 else
863 {
864 intExperIndex[ipLev] = -1;
865 }
866
867 }
868 else
869 {
870 //Stop looking for levels when the array is full
871 if( ipLev == nMolLevs )
872 break;
873
874 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
875 elvlcstream.get(obtemp,lvl_nrg_col);
876 if(atof(obtemp) != 0.)
877 {
878 fenergy = (realnum) atof(obtemp);
879 }
880 //else
881 //fprintf( ioQQQ," WARNING: Switched to theoretical energies for species %s, level %3li\n", dBaseSpecies[intNS].chLabel, ipLev );
882 dBaseStatesEnergy[ipLev] = fenergy;
883 }
884
885 elvlcstream.ignore(INT_MAX,'\n');
886
887 }
888 else
889 {
890 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename);
891 fclose( ioProtCollData );
893 }
894 }
895
896 for( long ipLev=1; ipLev<nMolLevs; ipLev++ )
897 {
898 /*To check if energy levels are in order*/
899 /*If the energy levels are not in order a flag is set*/
900 if (dBaseStatesEnergy[ipLev] < dBaseStatesEnergy[ipLev-1])
901 {
902 lgEneLevOrder = false;
903 }
904 }
905
906 // malloc vector for new energy-sorted indices
907 intNewIndex = (long int*)MALLOC((unsigned long)(nElvlcLines)* sizeof(long int));
908
909 if(lgEneLevOrder == false)
910 {
911 /*If the energy levels are not in order and the flag is set(lgEneLevOrder=FALSE)
912 then we sort the energy levels to find the relation matrix between the old and new indices*/
913 intOldIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int));
914 /*We now do a modified quick sort*/
915 spsort(dBaseStatesEnergy,nMolLevs,intOldIndex,2,&interror);
916 /*intNewIndex has the key to map*/
917 for( long i=0; i<nMolLevs; i++ )
918 {
919 ASSERT( intOldIndex[i] < nMolLevs );
920 intNewIndex[intOldIndex[i]] = i;
921 }
922
923 if( nMolLevs < nElvlcLines )
924 {
925 /* If any chianti energy levels do not have experimental values
926 * the vector that relates experimental levels to the default will not be the same size
927 * as the stored experimental energy levels.
928 * Therefore this code is required to pad the rest of the sorted vector.
929 * Any index larger than nMolLevs set to -1 */
930 for( long i=nMolLevs; i<nElvlcLines; i++)
931 {
932 intNewIndex[i] = -1;
933 }
934 }
935
936 if(DEBUGSTATE)
937 {
938 for( long i=0; i<nMolLevs; i++ )
939 {
940 printf("The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]);
941 }
942 for( long i=0; i<nMolLevs; i++ )
943 {
944 printf("The %ld value of the sorted energy array is %f \n",i,dBaseStatesEnergy[i]);
945 }
946 }
947 free( intOldIndex );
948 }
949 else
950 {
951 /* if energies already in correct order, new index is same as original. */
952 for( long i=0; i<nMolLevs; i++ )
953 {
954 intNewIndex[i] = i;
955 }
956
957 if( nMolLevs < nElvlcLines )
958 {
959 /*Pad the sorted vector so that it has the correct number of elements
960 * Any index larger than nMolLevs set to -1.*/
961 for( long i=nMolLevs; i<nElvlcLines; i++)
962 {
963 intNewIndex[i] = -1;
964 }
965 }
966 }
967
968 /* insist that there the intNewIndex values are unique */
969 for( long i=0; i<nMolLevs; i++ )
970 {
971 for( long j=i+1; j<nMolLevs; j++ )
972 {
973 ASSERT( intNewIndex[i] != intNewIndex[j] );
974 }
975 }
976
977 /*This will print out the relational matrix
978 * which includes the original index,the new experimental index,
979 * and the sorted experimental index.
980 * The numbers NOT be on the C scale, so 1 is the lowest level. */
981 if( DEBUGSTATE && atmdat.lgChiantiExp )
982 {
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++ )
986 {
987 if( intExperIndex[ipLevOld] != -1)
988 printf("%li\t%li\t%li\n",ipLevOld+1,intExperIndex[ipLevOld]+1,intNewIndex[intExperIndex[ipLevOld]]+1);
989 }
990 printf("\n");
991 }
992
993 //After we read in the energies we rewind the energy levels file
994 elvlcstream.seekg(0,ios::beg);
995
996 //lvl_skipto_statwt is the # of columns to skip to statwt from left
997 const int lvl_skipto_statwt = 37;
998 //lvl_statwt_col is the length (+1) of statwt
999 const int lvl_statwt_col = 4;
1000 //Read in stat weight and energy
1001 for( long ipLevOld=0; ipLevOld<nElvlcLines; ipLevOld++ )
1002 {
1003 long ipLevNew = 0;
1004 if( atmdat.lgChiantiExp )
1005 {
1006 /* We know that non-experimental levels are stored as -1
1007 * in the observed/experimental vector.
1008 * Use that to skip over those values. */
1009
1010 if( intExperIndex[ipLevOld] == -1 )
1011 {
1012 elvlcstream.ignore(INT_MAX,'\n');
1013 continue;
1014 }
1015 else
1016 {
1017 /* Store values based on the sorted experimental indices */
1018 ipLevNew = intNewIndex[intExperIndex[ipLevOld]];
1019 }
1020 }
1021 else
1022 {
1023 /* With level trimming on it is possible that there can be rows that
1024 * have to be skipped when using theoretical
1025 * since the levels no longer exist */
1026 if( intNewIndex[ipLevOld] == -1 )
1027 {
1028 elvlcstream.ignore(INT_MAX,'\n');
1029 continue;
1030 }
1031 else
1032 {
1033 ipLevNew = intNewIndex[ipLevOld];
1034 }
1035 }
1036
1037 char gtemp[lvl_statwt_col],thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
1038
1039 /*information needed for label*/
1040 strcpy(dBaseStates[intNS][ipLevNew].chLabel(), " ");
1041 strncpy(dBaseStates[intNS][ipLevNew].chLabel(),dBaseSpecies[intNS].chLabel, 4);
1042 dBaseStates[intNS][ipLevNew].chLabel()[4] = '\0';
1043
1044 //Move cursor to and read statwt
1045 elvlcstream.seekg(lvl_skipto_statwt,ios::cur);
1046 elvlcstream.get(gtemp,lvl_statwt_col);
1047 fstatwt = (realnum)atof(gtemp);
1048
1049 if(fstatwt <= 0.)
1050 {
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);
1054 fstatwt = 1.f;
1055 //cdEXIT(EXIT_FAILURE);
1056 }
1057 dBaseStates[intNS][ipLevNew].g() = fstatwt;
1058
1059 //Read nrg
1060 elvlcstream.get(thtemp,lvl_nrg_col);
1061
1062 fenergy = (realnum) atof(thtemp);
1063
1064 /* If we are looking for theoretical energies
1065 * move over a couple columns before reading in the energy*/
1066 if( !atmdat.lgChiantiExp )
1067 {
1068 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
1069 elvlcstream.get(obtemp,lvl_nrg_col);
1070 if(atof(obtemp) != 0.)
1071 {
1072 fenergy = (realnum) atof(obtemp);
1073 }
1074 }
1075 elvlcstream.ignore(INT_MAX,'\n');
1076 dBaseStates[intNS][ipLevNew].energy().set(fenergy,"cm^-1");
1077 }
1078 //Close the level stream
1079 elvlcstream.close();
1080
1081 // highest energy transition in chianti
1082 dBaseSpecies[intNS].maxWN = 0.;
1083 /* fill in all transition energies, can later overwrite for specific radiative transitions */
1084 for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
1085 tr!= dBaseTrans[intNS].end(); ++tr)
1086 {
1087 int ipHi = (*tr).ipHi();
1088 int ipLo = (*tr).ipLo();
1089 fenergyWN = (realnum)MAX2( ENERGY_MIN_WN , dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN() );
1090
1091 (*tr).EnergyWN() = fenergyWN;
1092 (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
1093 dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,fenergyWN);
1094 }
1095
1096 if(DEBUGSTATE)
1097 {
1098 /* Print out the stored data for each level.
1099 * Level index is NOT on C scale. */
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++ )
1103 {
1104 printf("%li\t%e\t%.1f\n",ipLo+1,dBaseStates[intNS][ipLo].energy().WN(),dBaseStates[intNS][ipLo].g());
1105 }
1106 }
1107
1108 /************************************************************************/
1109 /*** Read in the level numbers, Einstein A and transition wavelength ***/
1110 /************************************************************************/
1111
1112 //Count the number of rows first
1113 long wgfarows = -1;
1114 if (wgfastream.is_open())
1115 {
1116 int nj = 0;
1117 char otemp[eof_col];
1118 while(nj != -1)
1119 {
1120 wgfastream.get(otemp,eof_col);
1121 wgfastream.ignore(INT_MAX,'\n');
1122 nj = atoi(otemp);
1123 wgfarows++;
1124 }
1125 wgfastream.seekg(0,ios::beg);
1126 }
1127 else
1128 fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
1129
1130 //line_index_col is the length(+1) of the level indexes in the WGFA file
1131 const int line_index_col = 6;
1132 //line_nrg_to_eina is the # of columns to skip from wavelength to eina in WGFA file
1133 const int line_nrg_to_eina =15;
1134 //line_eina_col is the length(+1) of einsteinA in WGFA
1135 const int line_eina_col = 16;
1136 char lvltemp[line_index_col];
1137 //Start reading WGFA file
1138 if (wgfastream.is_open())
1139 {
1140 for (long i = 0;i<wgfarows;i++)
1141 {
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;
1146
1147 if( atmdat.lgChiantiExp )
1148 {
1149 /* If either upper or lower index is -1 in the relational vector,
1150 * skip that line in the wgfa file.
1151 * Otherwise translate the level indices.*/
1152 if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1153 {
1154 wgfastream.ignore(INT_MAX,'\n');
1155 continue;
1156 }
1157 else
1158 {
1159 ipLo = intNewIndex[intExperIndex[ipLo]];
1160 ipHi = intNewIndex[intExperIndex[ipHi]];
1161 }
1162 }
1163 else
1164 {
1165 /* With level trimming on it is possible that there can be rows that
1166 * have to be skipped when using theoretical
1167 * since the levels no longer exist */
1168 if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
1169 {
1170 wgfastream.ignore(INT_MAX,'\n');
1171 continue;
1172 }
1173 else
1174 {
1175 ipLo = intNewIndex[ipLo];
1176 ipHi = intNewIndex[ipHi];
1177 }
1178 }
1179
1180 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1181 {
1182 // skip these lines
1183 wgfastream.ignore(INT_MAX,'\n');
1184 continue;
1185 }
1186
1187 if( ipHi == ipLo )
1188 {
1189 fprintf( ioQQQ," WARNING: Upper level = lower for a radiative transition in %s. Ignoring.\n", chTraFilename );
1190 wgfastream.ignore(INT_MAX,'\n');
1191 continue;
1192 }
1193
1194
1195 ASSERT( ipHi != ipLo );
1196 ASSERT(ipHi >= 0);
1197 ASSERT(ipLo >= 0);
1198
1199 // sometimes the CHIANTI datafiles list the highest index first as in the middle of these five lines in ne_10.wgfa:
1200 // ...
1201 // 8 10 187.5299 0.000e+00 4.127e+05 3d 2D1.5 - 4s 2S0.5 E2
1202 // 9 10 187.6573 0.000e+00 6.197e+05 3d 2D2.5 - 4s 2S0.5 E2
1203 // 11 10 4842624.0000 1.499e-05 9.423e-06 4p 2P0.5 - 4s 2S0.5 E1
1204 // 1 11 9.7085 1.892e-02 6.695e+11 1s 2S0.5 - 4p 2P0.5 E1
1205 // 2 11 48.5157 6.787e-02 9.618e+10 2s 2S0.5 - 4p 2P0.5 E1
1206 // ...
1207 // so, just set ipHi (ipLo) equal to the max (min) of the two indices.
1208 // NB NB NB it looks like this may depend upon whether one uses observed or theoretical energies.
1209
1210 //Read in wavelengths
1211 char trantemp[lvl_nrg_col];
1212 wgfastream.get(trantemp,lvl_nrg_col);
1213 fWLAng = (realnum)atof(trantemp);
1214
1215 /* \todo 2 CHIANTI labels the H 1 2-photon transition as z wavelength of zero.
1216 * Should we just ignore all of the wavelengths in this file and use the
1217 * difference of level energies instead. */
1218
1219 if( atmdat.lgChiantiExp )
1220 {
1221 /* Sometimes the wgfa file lists the levels to where ipLo > ipHi.
1222 * This seems to be related to the theoretical energies being out of order for those transitions.
1223 * When this is true it seems that the associated wavelength is negative which means they are using theoretical energies,
1224 * even though they have experimental energies.
1225 * For now, print that the levels are switched and ignore the lines if the wavelength is negative. */
1226
1227 if( ipHi < ipLo )
1228 {
1229 if( strcmp(dBaseSpecies[intNS].chLabel,"Fe 3") == 0)
1230 {
1231 long ipLoTemp = ipLo;
1232 long ipHiTemp = ipHi;
1233 ipHi = max( ipHiTemp, ipLoTemp );
1234 ipLo = min( ipHiTemp, ipLoTemp );
1235 if( atmdat.lgChiantiPrint )
1236 {
1237 fprintf( ioQQQ," WARNING: Swapped level indices for species %6s, indices %3li %3li with wavelength %e \n",
1238 dBaseSpecies[intNS].chLabel, ipLoTemp, ipHiTemp,fWLAng );
1239 }
1240 }
1241 else if( atmdat.lgChiantiPrint )
1242 {
1243 fprintf( ioQQQ," WARNING: Upper and Lower Indices are reversed.Species: %6s, Indices: %3li %3li Wavelength: %e \n",
1244 dBaseSpecies[intNS].chLabel, ipLo+1, ipHi+1,fWLAng );
1245 }
1246 }
1247
1248 if( fWLAng < 0.)
1249 {
1250 // skip these lines
1251 wgfastream.ignore(INT_MAX,'\n');
1252 if( atmdat.lgChiantiPrint )
1253 {
1254 fprintf(ioQQQ,"WARNING: Skipping Transition %li to %li of %s.\n",ipLo+1,ipHi+1,dBaseSpecies[intNS].chLabel);
1255 }
1256 continue;
1257 }
1258
1259 }
1260 else
1261 {
1262 /* If Cloudy is using theoretical energies, then just flip the levels. */
1263 if( ipHi < ipLo )
1264 {
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 );
1271 }
1272 }
1273
1274 /* If the given wavelength is negative, then theoretical enegies are being used.
1275 * Take the difference in stored theoretical energies.
1276 * It should equal the absolute value of the wavelength in the wgfa file. */
1277 if( !atmdat.lgChiantiExp && fWLAng <= 0. )
1278 {
1279 //if( fWLAng < 0.)
1280 //fprintf( ioQQQ," WARNING: Negative wavelength for species %6s, indices %3li %3li \n", dBaseSpecies[intNS].chLabel, ipLo, ipHi);
1281 fWLAng = (realnum)(1e8/
1282 (dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN()));
1283 }
1284 //Skip from end of Wavelength to Einstein A and read in
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 )
1289 {
1290 static bool notPrintedYet = true;
1291 if( notPrintedYet && atmdat.lgChiantiPrint)
1292 {
1293 fprintf( ioQQQ," WARNING: Radiative rate(s) equal to zero in %s.\n", chTraFilename );
1294 notPrintedYet = false;
1295 }
1296 wgfastream.ignore(INT_MAX,'\n');
1297 continue;
1298 }
1299
1300 fixit(); // may need to do something with these rates
1301 //Read in the rest of the line and look for auto
1302 wgfastream.getline(chLine,INT_MAX);
1303 TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
1304 if( nMatch("auto", chLine) )
1305 {
1306 if( (*tr).hasEmis() )
1307 {
1308 (*tr).Emis().AutoIonizFrac() =
1309 feinsteina/((*tr).Emis().Aul() + feinsteina);
1310 ASSERT( (*tr).Emis().AutoIonizFrac() >= 0. );
1311 ASSERT( (*tr).Emis().AutoIonizFrac() <= 1. );
1312 }
1313 continue;
1314 }
1315
1316 if( (*tr).hasEmis() )
1317 {
1318 fprintf(ioQQQ," PROBLEM duplicate transition found by atmdat_chianti in %s, "
1319 "wavelength=%f\n", chTraFilename,fWLAng);
1320 fclose( ioProtCollData );
1322 }
1323 (*tr).AddLine2Stack();
1324 (*tr).Emis().Aul() = feinsteina;
1325 (*tr).WLAng() = fWLAng;
1326
1327 fenergyWN = (realnum)(1e+8/fWLAng);
1328
1329 // TODO::Check the wavelength in the file with the difference in energy levels
1330
1331 (*tr).EnergyWN() = fenergyWN;
1332 (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
1333 (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
1334
1335 if(DEBUGSTATE)
1336 {
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() );
1341 }
1342 }
1343 }
1344 else fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
1345 wgfastream.close();
1346
1347 /* Malloc space for splines */
1348 AtmolCollSplines[intNS] = (CollSplinesArray***)MALLOC(nMolLevs *sizeof(CollSplinesArray**));
1349 for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
1350 {
1351 AtmolCollSplines[intNS][ipHi] = (CollSplinesArray **)MALLOC((unsigned long)(nMolLevs)*sizeof(CollSplinesArray *));
1352 for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
1353 {
1354 AtmolCollSplines[intNS][ipHi][ipLo] =
1355 (CollSplinesArray *)MALLOC((unsigned long)(ipNCOLLIDER)*sizeof(CollSplinesArray ));
1356
1357 for( long k=0; k<ipNCOLLIDER; k++ )
1358 {
1359 /* initialize all spline variables */
1360 AtmolCollSplines[intNS][ipHi][ipLo][k].collspline = NULL;
1361 AtmolCollSplines[intNS][ipHi][ipLo][k].SplineSecDer = NULL;
1362 AtmolCollSplines[intNS][ipHi][ipLo][k].nSplinePts = -1;
1363 AtmolCollSplines[intNS][ipHi][ipLo][k].intTranType = -1;
1364 AtmolCollSplines[intNS][ipHi][ipLo][k].EnergyDiff = BIGDOUBLE;
1365 AtmolCollSplines[intNS][ipHi][ipLo][k].ScalingParam = BIGDOUBLE;
1366 }
1367 }
1368 }
1369
1370 /************************************/
1371 /*** Read in the collisional data ***/
1372 /************************************/
1373
1374 // ipCollider 0 is electrons
1375 for( long ipCollider=0; ipCollider < 1; ipCollider++ )
1376 {
1377 char chFilename[FILENAME_PATH_LENGTH_2];
1378
1379 if( ipCollider == ipELECTRON )
1380 {
1381 strcpy( chFilename, chEleColFilename );
1382 }
1383 else if( ipCollider == ipPROTON )
1384 {
1385 if( !lgProtonData )
1386 break;
1387 fprintf( ioQQQ," WARNING: Chianti proton collision data have different format than electron data files!\n" );
1388 strcpy( chFilename, chProColFilename );
1389 }
1390 else
1391 TotalInsanity();
1392
1393 /*Dummy string used for convenience*/
1394 strcpy(chLine,"A");
1395
1396 fstream splupsstream;
1397 open_data(splupsstream, chFilename,mode_r);
1398
1399 //cs_eof_col is the length(+1) of the first column used for finding the end of file
1400 const int cs_eof_col = 4;
1401 //cs_index_col is the length(+1) of the indexes in the CS file
1402 const int cs_index_col = 4;
1403 //cs_trantype_col is the length(+1) of the transition type in the CS file
1404 const int cs_trantype_col = 4;
1405 //cs_values_col is the length(+1) of the other values in the CS file
1406 //including: GF, nrg diff, scaling parameter, and spline points
1407 const int cs_values_col = 11;
1408 //Determine the number of rows in the CS file
1409 if (splupsstream.is_open())
1410 {
1411 int nj = 0;
1412 //splupslines is -1 since the loop runs 1 extra time
1413 long splupslines = -1;
1414 char otemp[cs_eof_col];
1415 while(nj != -1)
1416 {
1417 splupsstream.get(otemp,cs_eof_col);
1418 splupsstream.ignore(INT_MAX,'\n');
1419 nj = atoi(otemp);
1420 splupslines++;
1421 }
1422 splupsstream.seekg(0,ios::beg);
1423
1424 for (int m = 0;m<splupslines;m++)
1425 {
1426 if( ipCollider == ipELECTRON )
1427 {
1428 splupsstream.seekg(6,ios::cur);
1429 }
1430
1431 /* level indices */
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;
1436
1437 /* If either upper or lower index is -1 in the relational vector,
1438 * skip that line in the splups file.
1439 * Otherwise translate the level indices.*/
1440 if( atmdat.lgChiantiExp )
1441 {
1442 if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1443 {
1444 splupsstream.ignore(INT_MAX,'\n');
1445 continue;
1446 }
1447 else
1448 {
1449 ipLo = intNewIndex[intExperIndex[ipLo]];
1450 ipHi = intNewIndex[intExperIndex[ipHi]];
1451 }
1452
1453 }
1454 else
1455 {
1456 /* With level trimming on it is possible that there can be rows that
1457 * have to be skipped when using theoretical
1458 * since the levels no longer exist */
1459 if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
1460 {
1461 splupsstream.ignore(INT_MAX,'\n');
1462 continue;
1463 }
1464 else
1465 {
1466 ipLo = intNewIndex[ipLo];
1467 ipHi = intNewIndex[ipHi];
1468 }
1469 }
1470
1471 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1472 {
1473 // skip these transitions
1474 splupsstream.ignore(INT_MAX,'\n');
1475 continue;
1476 }
1477
1478 ASSERT( ipLo < nMolLevs );
1479 ASSERT( ipHi < nMolLevs );
1480 /*Transition Type*/
1481 splupsstream.get(otemp,cs_trantype_col);
1482 intTranType = atoi(otemp);
1483 char qtemp[cs_values_col];
1484 splupsstream.get(qtemp,cs_values_col);
1485 /*Energy Difference*/
1486 splupsstream.get(qtemp,cs_values_col);
1487 fEnergyDiff = atof(qtemp);
1488 /*Scaling Parameter*/
1489 splupsstream.get(qtemp,cs_values_col);
1490 fScalingParam = atof(qtemp);
1491
1492 ASSERT( ipLo < nMolLevs );
1493 ASSERT( ipHi < nMolLevs );
1494 ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline == NULL );
1495 ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer == NULL );
1496
1497 const int CHIANTI_SPLINE_MAX=9, CHIANTI_SPLINE_MIN=5;
1498 STATIC_ASSERT(CHIANTI_SPLINE_MAX > CHIANTI_SPLINE_MIN);
1499
1500 /*We malloc the space here*/
1501 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline =
1502 (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
1503 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer =
1504 (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
1505
1506 /* always read at least CHIANTI_SPLINE_MIN */
1507 for( intsplinepts=0; intsplinepts<=CHIANTI_SPLINE_MAX; intsplinepts++ )
1508 {
1509 //Look at the next character to see if it is the end of line.
1510 char p = splupsstream.peek();
1511 if( p == '\n' )
1512 {
1513 /* there should be 5 or 9 spline points. If we got EOL,
1514 * insist that we were trying to read the 6th or 10th. */
1515 if( (intsplinepts != CHIANTI_SPLINE_MAX) && (intsplinepts != CHIANTI_SPLINE_MIN) )
1516 {
1517 static bool notPrintedYet = true;
1518 if( notPrintedYet )
1519 {
1520 fprintf( ioQQQ, " WARNING: Wrong number of spline points in %s.\n", chFilename);
1521 notPrintedYet = false;
1522 }
1523 for( long j=intsplinepts-1; j<CHIANTI_SPLINE_MAX; j++ )
1524 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[j] = 0.;
1525 }
1526 else
1527 break;
1528 }
1529 else
1530 {
1531 if( intsplinepts >= CHIANTI_SPLINE_MAX )
1532 {
1533 fprintf( ioQQQ, " WARNING: More spline points than expected in %s, indices %3li %3li. Ignoring extras.\n", chFilename, ipHi, ipLo );
1534 break;
1535 }
1536 ASSERT( intsplinepts < CHIANTI_SPLINE_MAX );
1537 double temp;
1538 //Store a single spline point then look for more
1539 splupsstream.get(qtemp,cs_values_col);
1540 temp = atof(qtemp);
1541 if(temp < 0)
1542 {
1543 temp = 0.;
1544 }
1545 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intsplinepts] = temp;
1546 }
1547
1548 }
1549
1550 ASSERT( (intsplinepts == CHIANTI_SPLINE_MAX) || (intsplinepts == CHIANTI_SPLINE_MIN) );
1551
1552 /*The zeroth element contains the number of spline points*/
1553 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].nSplinePts = intsplinepts;
1554 /*Transition type*/
1555 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].intTranType = intTranType;
1556 /*Energy difference between two levels in Rydbergs*/
1557 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].EnergyDiff = fEnergyDiff;
1558 /*Scaling parameter C*/
1559 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].ScalingParam = fScalingParam;
1560
1561 /*Once the spline points have been filled,fill the second derivatives structure*/
1562 /*Creating spline points array*/
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));
1566
1567 // should be able to just loop to intsplinepts -- ASSERT above protects
1568 if(intsplinepts == CHIANTI_SPLINE_MIN)
1569 {
1570 for(intxs=0;intxs<CHIANTI_SPLINE_MIN;intxs++)
1571 {
1572 xs[intxs] = 0.25*intxs;
1573 spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
1574 }
1575 }
1576 else if(intsplinepts == CHIANTI_SPLINE_MAX)
1577 {
1578 for(intxs=0;intxs<CHIANTI_SPLINE_MAX;intxs++)
1579 {
1580 xs[intxs] = 0.125*intxs;
1581 spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
1582 }
1583 }
1584 else
1585 TotalInsanity();
1586
1587 spline(xs, spl,intsplinepts,2e31,2e31,spl2);
1588
1589 /*Filling the second derivative structure*/
1590 for( long i=0; i<intsplinepts; i++)
1591 {
1592 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer[i] = spl2[i];
1593 }
1594
1595 free( xs );
1596 free( spl );
1597 free( spl2 );
1598 splupsstream.ignore(INT_MAX,'\n');
1599 }
1600 splupsstream.close();
1601 }
1602 }
1603
1604 // free some memory
1605 free( dBaseStatesEnergy );
1606 free( intNewIndex );
1607 if( atmdat.lgChiantiExp)
1608 {
1609 free( intExperIndex );
1610 }
1611
1612 // close open file handles
1613 fclose( ioElecCollData );
1614 if( lgProtonData )
1615 fclose( ioProtCollData );
1616
1617 return;
1618}
t_atmdat atmdat
Definition atmdat.cpp:6
const double ENERGY_MIN_WN
void atmdat_CHIANTI_readin(long intNS, char *chPrefix)
void atmdat_STOUT_readin(long intNS, char *chPrefix)
#define DEBUGSTATE
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
struct t_CollSplinesArray CollSplinesArray
long nMatch(const char *chKey, const char *chCard)
Definition service.cpp:451
char tolower(char c)
Definition cddefines.h:691
#define MALLOC(exp)
Definition cddefines.h:501
#define MIN3(a, b, c)
Definition cddefines.h:766
struct t_StoutColls StoutColls
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long max(int a, long b)
Definition cddefines.h:775
#define STATIC_ASSERT(x)
Definition cddefines.h:113
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition service.cpp:1100
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
void uncaps(char *chCard)
Definition service.cpp:263
void fixit(void)
Definition service.cpp:991
TransitionProxy::iterator iterator
Definition transition.h:280
@ ipPROTON
Definition collision.h:10
@ ipNCOLLIDER
Definition collision.h:18
@ ipELECTRON
Definition collision.h:9
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const ios_base::openmode mode_r
Definition cpu.h:212
const double BIGDOUBLE
Definition cpu.h:194
double RefIndex(double EnergyWN)
double GetGF(double trans_prob, double enercm, double gup)
static double * g
Definition species2.cpp:28
vector< qList > dBaseStates
Definition taulines.cpp:15
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
vector< multi_arr< int, 2 > > ipdBaseTrans
Definition taulines.cpp:16
StoutColls **** StoutCollData
Definition taulines.cpp:20
species * dBaseSpecies
Definition taulines.cpp:14
CollSplinesArray **** AtmolCollSplines
Definition taulines.cpp:19
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
Definition thirdparty.h:117
t_trace trace
Definition trace.cpp:5