cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_chianti.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 
11 #define DEBUGSTATE false
12 
13 void atmdat_Chianti_readin( long intNS )
14 {
15  DEBUG_ENTRY( "atmdat_Chianti_readin()" );
16 
17  int intCS,intCollIndex = -10000,intsplinepts,intTranType,intxs;
18  /* type is set to 0 for non chianti and 1 for chianti*/
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;
23 
24  //Energy Column Start
25  //Statistical Weight Start
26  //Lower Level Start
27  //Upper Level Start
28  //Einstein A Start
29  //Wavelength of transition start
30  long ECS,SWS,LLS,ULS,EAS,WOTS;
31 
32  char chLine[FILENAME_PATH_LENGTH_2] ,
33  chEnFilename[FILENAME_PATH_LENGTH_2],
34  chTraFilename[FILENAME_PATH_LENGTH_2],
35  chEleColFilename[FILENAME_PATH_LENGTH_2],
36  chProColFilename[FILENAME_PATH_LENGTH_2];
37 
38  char *chDesired;
39  realnum *atmolStatesEnergy;
40  long int *intNewIndex=NULL,*intOldIndex=NULL;
41  int interror;
42  bool lgProtonData=false,lgEneLevOrder;
43 
44  /*This is to identify where the fields start */
45  /*In Chianti(1),the energy field in cm-1 starts at 40 and has a width 15*/
46  ECS = 71;/*We use the theoretical energy in cm-1*/
47  SWS = 38;
48  LLS = 1;
49  ULS = 6;
50  EAS = 41;
51  WOTS = 11;
52 
53  /*For the CHIANTI DATABASE*/
54  /*First Opening the energy levels file*/
55  strcpy( chEnFilename , Species[intNS].chptrSpName );
56  strcat( chEnFilename , ".elvlc");
57  uncaps( chEnFilename );
58  if(DEBUGSTATE)
59  printf("The data file is %s \n",chEnFilename);
60 
61  /*Open the files*/
62  if( trace.lgTrace )
63  fprintf( ioQQQ," moldat_readin opening %s:",chEnFilename);
64 
65  atmolLevDATA = open_data( chEnFilename, "r" );
66 
67  /*Second:Opening the transition probabilities file*/
68  strcpy( chTraFilename , Species[intNS].chptrSpName );
69  strcat( chTraFilename , ".wgfa");
70  uncaps( chTraFilename );
71  if(DEBUGSTATE)
72  printf("The data file is %s \n",chTraFilename);
73 
74  /*Open the files*/
75  if( trace.lgTrace )
76  fprintf( ioQQQ," moldat_readin opening %s:",chTraFilename);
77 
78  atmolTraDATA = open_data( chTraFilename, "r" );
79 
80  /*Third: Opening the electron collision data*/
81  strcpy( chEleColFilename , Species[intNS].chptrSpName );
82  strcat( chEleColFilename , ".splups");
83  uncaps( chEleColFilename );
84  if(DEBUGSTATE)
85  printf("The data file is %s \n",chEleColFilename);
86  /*Open the files*/
87  if( trace.lgTrace )
88  fprintf( ioQQQ," moldat_readin opening %s:",chEleColFilename);
89 
90  atmolEleColDATA = open_data( chEleColFilename, "r" );
91 
92  /*Fourth: Opening the proton collision data*/
93  strcpy( chProColFilename , Species[intNS].chptrSpName );
94  strcat( chProColFilename , ".psplups");
95  uncaps( chProColFilename );
96  if(DEBUGSTATE)
97  printf("The data file is %s \n",chProColFilename);
98  /*Open the files*/
99  if( trace.lgTrace )
100  fprintf( ioQQQ," moldat_readin opening %s:",chProColFilename);
101  /*We will set a flag here to indicate if the proton collision strengths are available */
102  if( ( atmolProColDATA = fopen( chProColFilename , "r" ) ) != NULL )
103  {
104  lgProtonData = true;
105  fclose( atmolProColDATA );
106  atmolProColDATA = NULL;
107  }
108  else
109  {
110  lgProtonData = false;
111  }
112 
113  nMolLevs = 0;
114  lgEneLevOrder = true;
115  while( read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) != NULL )
116  {
117  chDesired = strtok(chLine," \t");
118  intCS = atoi(chDesired);
119  if (intCS == -1)
120  {
121  break;
122  }
123  else
124  nMolLevs++;
125  }
126 
127  if(DEBUGSTATE)
128  printf("The number of energy levels is %li",nMolLevs);
129 
130  if( nMolLevs <= 0 )
131  {
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" );
134  cdEXIT( EXIT_FAILURE );
135  }
136 
137  Species[intNS].numLevels_max = nMolLevs;
138  Species[intNS].numLevels_local = nMolLevs;
139  /*Malloc the energy array to check if the energy levels are in order*/
140  atmolStatesEnergy = (realnum*)MALLOC((unsigned long)(nMolLevs)*sizeof(realnum));
141  /*malloc the States array*/
142  atmolStates[intNS] = (quantumState *)MALLOC((size_t)(nMolLevs)*sizeof(quantumState));
143  /*malloc the Transition array*/
144  atmolTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)nMolLevs);
145  atmolTrans[intNS][0] = NULL;
146  for( ipHi = 1; ipHi < nMolLevs; ipHi++)
147  {
148  atmolTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi);
149  for(ipLo = 0; ipLo < ipHi; ipLo++)
150  {
151  atmolTrans[intNS][ipHi][ipLo].Lo = &atmolStates[intNS][ipLo];
152  atmolTrans[intNS][ipHi][ipLo].Hi = &atmolStates[intNS][ipHi];
153  atmolTrans[intNS][ipHi][ipLo].Emis = &DummyEmis;
154  }
155  }
156 
157  for( intcollindex = 0; intcollindex<NUM_COLLIDERS; intcollindex++ )
158  {
159  CollRatesArray[intNS][intcollindex] = NULL;
160  }
161  /*Allocating space just for the electron*/
162  CollRatesArray[intNS][0] = (double**)MALLOC((nMolLevs)*sizeof(double*));
163  for( ipHi = 0; ipHi<nMolLevs; ipHi++ )
164  {
165  CollRatesArray[intNS][0][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double));
166  for( ipLo = 0; ipLo<nMolLevs; ipLo++)
167  {
168  CollRatesArray[intNS][0][ipHi][ipLo] = 0.0;
169  }
170  }
171  /*Allocating space for the proton*/
172  if(lgProtonData)
173  {
174  CollRatesArray[intNS][1] = (double**)MALLOC((nMolLevs)*sizeof(double*));
175  for( ipHi = 0; ipHi<nMolLevs; ipHi++ )
176  {
177  CollRatesArray[intNS][1][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double));
178  for( ipLo = 0; ipLo<nMolLevs; ipLo++)
179  {
180  CollRatesArray[intNS][1][ipHi][ipLo] = 0.0;
181  }
182  }
183  }
184  /*Rewind the energy levels files*/
185  if( fseek( atmolLevDATA , 0 , SEEK_SET ) != 0 )
186  {
187  fprintf( ioQQQ, " moldat_readin could not rewind %s.\n",chEnFilename);
188  cdEXIT(EXIT_FAILURE);
189  }
190 
191 
192  /*Reading in the energy and checking that they are in order*/
193  for( i=0; i<nMolLevs; i++ )
194  {
195  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
196  {
197  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename);
198  cdEXIT(EXIT_FAILURE);
199  }
200  fenergy = (realnum) atof(&chLine[ECS]);
201  atmolStatesEnergy[i] = fenergy;
202 
203  /*To check if energy levels are in order*/
204  /*If the energy levels are not in order a flag is set*/
205  if(i>0)
206  {
207  if (atmolStatesEnergy[i] < atmolStatesEnergy[i-1])
208  {
209  lgEneLevOrder = false;
210  }
211  }
212  }
213 
214  if(DEBUGSTATE)
215  {
216  for(i=0;i<nMolLevs;i++)
217  {
218  printf("The %ld value of the unsorted energy array is %f \n",i,atmolStatesEnergy[i]);
219  }
220  }
221 
222  intNewIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int));
223 
224  if(lgEneLevOrder == false)
225  {
226  /*If the energy levels are not in order and the flag is set(lgEneLevOrder=FALSE)
227  then we sort the energy levels to find the relation matrix between the old and new indices*/
228  intOldIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int));
229  /*We now do a modified quick sort*/
230  spsort(atmolStatesEnergy,nMolLevs,intOldIndex,2,&interror);
231  /*intNewIndex has the key to map*/
232  for( i=0; i<nMolLevs; i++ )
233  {
234  ASSERT( intOldIndex[i] < nMolLevs );
235  intNewIndex[intOldIndex[i]] = i;
236  }
237  if(DEBUGSTATE)
238  {
239  for(i=0;i<nMolLevs;i++)
240  {
241  printf("The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]);
242  }
243  for( i=0; i<nMolLevs; i++ )
244  {
245  printf("The %ld value of the sorted energy array is %f \n",i,atmolStatesEnergy[i]);
246  }
247  }
248 
249  free( intOldIndex );
250  }
251  else
252  {
253  /* if energies already in correct order, new index is same as original. */
254  for( i=0; i<nMolLevs; i++ )
255  {
256  intNewIndex[i] = i;
257  }
258  }
259 
260  /* insist that there the intNewIndex values are unique */
261  for( i=0; i<nMolLevs; i++ )
262  {
263  for( j=i+1; j<nMolLevs; j++ )
264  {
265  ASSERT( intNewIndex[i] != intNewIndex[j] );
266  }
267  }
268 
269 
270  /*After we read in the energies we rewind the energy levels file */
271  if( fseek( atmolLevDATA , 0 , SEEK_SET ) != 0 )
272  {
273  fprintf( ioQQQ, " moldat_readin could not rewind %s.\n",chEnFilename);
274  cdEXIT(EXIT_FAILURE);
275  }
276 
277  for( i=0; i<nMolLevs; i++ )
278  {
279  long j = intNewIndex[i];
280 
281  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
282  {
283  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename);
284  cdEXIT(EXIT_FAILURE);
285  }
286  /*information needed for label*/
287  strcpy(atmolStates[intNS][j].chLabel, " ");
288  strncpy(atmolStates[intNS][j].chLabel,Species[intNS].chptrSpName, 4);
289  atmolStates[intNS][j].chLabel[4] = '\0';
290 
291  fstatwt = (realnum)atof(&chLine[SWS]);
292  if(fstatwt <= 0.)
293  {
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);
296  cdEXIT(EXIT_FAILURE);
297 
298  }
299  atmolStates[intNS][j].g = fstatwt;
300  fenergy = (realnum) atof(&chLine[ECS]);
301  atmolStates[intNS][j].energy = fenergy;
302 
303  if(DEBUGSTATE)
304  {
305  fprintf( ioQQQ, "The %ld value of g after populating is %f \n",
306  j,atmolStates[intNS][j].g);
307  fprintf( ioQQQ, "The %ld value of energy after populating is %f \n",
308  j,atmolStates[intNS][j].energy);
309  }
310  }
311 
312 
313  /* fill in all transition energies, can later overwrite for specific radiative transitions */
314  for( i=1; i<nMolLevs; i++ )
315  {
316  for( j=0; j<i; j++ )
317  {
318  fenergyWN = MAX2( (realnum)1E-10, atmolStates[intNS][i].energy - atmolStates[intNS][j].energy );
319  fenergyK = (realnum)(fenergyWN*T1CM);
320 
321  atmolTrans[intNS][i][j].EnergyK = fenergyK;
322  atmolTrans[intNS][i][j].EnergyWN = fenergyWN;
323  atmolTrans[intNS][i][j].EnergyErg = (realnum)ERG1CM *fenergyWN;
324  atmolTrans[intNS][i][j].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
325  }
326  }
327 
328  /************************************************************************/
329  /*** Read in the level numbers, Einstein A and transition wavelength ***/
330  /************************************************************************/
331  if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolTraDATA ) == NULL )
332  {
333  fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
334  cdEXIT(EXIT_FAILURE);
335  }
336 
337  while(atoi(chLine)!= -1)
338  {
339  ipLo = atoi(&chLine[LLS])-1;
340  ipHi = atoi(&chLine[ULS])-1;
341 
342  /* translate to properly sorted indices */
343  ipLo = intNewIndex[ipLo];
344  ipHi = intNewIndex[ipHi];
345 
346  ASSERT( ipHi > ipLo );
347 
348  fWLAng = (realnum)atof(&chLine[WOTS]);
349 
350  /* \todo 2 Chianti labels the H 1 2-photon transition as z wavelength of zero.
351  * Should we just ignore all of the wavelengths in this file and use the
352  * difference of level energies instead. */
353  if( fWLAng == 0. )
354  {
355  fWLAng = (realnum)1e8/
356  (atmolStates[intNS][ipHi].energy - atmolStates[intNS][ipLo].energy);
357  }
358 
359  feinsteina = (realnum)atof(&chLine[EAS]);
360  atmolTrans[intNS][ipHi][ipLo].Emis
361  = AddLine2Stack(feinsteina, &atmolTrans[intNS][ipHi][ipLo]);
362  atmolTrans[intNS][ipHi][ipLo].WLAng = fWLAng;
363  fenergyWN = (realnum)(1e+8/fWLAng);
364  fenergyK = (realnum)(fenergyWN*T1CM);
365 
366  atmolTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
367  atmolTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
368  atmolTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
369  atmolTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
370 
371  ASSERT( !isnan( atmolTrans[intNS][ipHi][ipLo].EnergyK ) );
372 
373  if(DEBUGSTATE)
374  {
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 );
379  }
380 
381  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolTraDATA ) == NULL )
382  {
383  fprintf( ioQQQ, " Data is expected on this line.\n");
384  fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
385  cdEXIT(EXIT_FAILURE);
386  }
387  }
388 
389  /* Malloc space for splines */
390 
391  AtmolCollSplines[intNS] = (CollSplinesArray***)MALLOC(nMolLevs *sizeof(CollSplinesArray**));
392  for( i=0; i<nMolLevs; i++ )
393  {
394  AtmolCollSplines[intNS][i] = (CollSplinesArray **)MALLOC((unsigned long)(nMolLevs)*sizeof(CollSplinesArray *));
395  for( j=0; j<nMolLevs; j++ )
396  {
397  AtmolCollSplines[intNS][i][j] =
398  (CollSplinesArray *)MALLOC((unsigned long)(NUM_COLLIDERS)*sizeof(CollSplinesArray ));
399 
400  for( long k=0; k<NUM_COLLIDERS; k++ )
401  {
402  /* initialize all spline variables */
403  AtmolCollSplines[intNS][i][j][k].collspline = NULL;
404  AtmolCollSplines[intNS][i][j][k].SplineSecDer = NULL;
405  AtmolCollSplines[intNS][i][j][k].nSplinePts = -1;
406  AtmolCollSplines[intNS][i][j][k].intTranType = -1;
407  AtmolCollSplines[intNS][i][j][k].gf = BIGDOUBLE;
408  AtmolCollSplines[intNS][i][j][k].EnergyDiff = BIGDOUBLE;
409  AtmolCollSplines[intNS][i][j][k].ScalingParam = BIGDOUBLE;
410  }
411  }
412  }
413 
414  fixit(); /* This is where the loop of colliders should begin */
415 
416  /*********************************************/
417  /*** Read in the electron collisional data ***/
418  /*********************************************/
419  intCollIndex = 0;
420  intCollTran = 0;
421  if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolEleColDATA ) == NULL )
422  {
423  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename);
424  cdEXIT(EXIT_FAILURE);
425  }
426  /*We need to find the number of collisional transitions */
427  while(atoi(chLine)!= -1)
428  {
429  intCollTran++;
430  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolEleColDATA ) == NULL )
431  {
432  fprintf( ioQQQ, " Data is expected on this line.\n");
433  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename);
434  cdEXIT(EXIT_FAILURE);
435 
436  }
437  }
438 
439  if(DEBUGSTATE)
440  printf("\n The number of collisional transitions is %li",intCollTran);
441 
442  /*We rewind the file*/
443  if( fseek( atmolEleColDATA , 0 , SEEK_SET ) != 0 )
444  {
445  fprintf( ioQQQ, " moldat_readin could not rewind %s.\n",chEleColFilename);
446  cdEXIT(EXIT_FAILURE);
447  }
448  /*Dummy string used for convenience*/
449  strcpy(chLine,"A");
450 
451  if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolEleColDATA ) == NULL )
452  {
453  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename);
454  cdEXIT(EXIT_FAILURE);
455  }
456 
457  /*If the file exists,the first line cannot be "-1"*/
458  if(atoi(chLine)== -1)
459  {
460  fprintf( ioQQQ, " If the file %s exists,the first line cannot be -1 .\n",chEleColFilename);
461  cdEXIT(EXIT_FAILURE);
462  }
463 
464  while(atoi(chLine)!= -1)
465  {
466  i = 1;
467  bool lgEOL;
468 
469  (void)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); // intAtomicNo
470  (void)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); // intIonStage
471  /* level indices */
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 );
475  /*gf*/
476  fGF = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
477  /*Energy Difference*/
478  fEnergyDiff = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
479  /*Scaling Parameter*/
480  fScalingParam = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
481 
482  /* translate to properly sorted indices */
483  ipLo = intNewIndex[ipLo];
484  ipHi = intNewIndex[ipHi];
485 
486  ASSERT( ipLo < nMolLevs );
487  ASSERT( ipHi < nMolLevs );
488  ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline == NULL );
489  ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].SplineSecDer == NULL );
490 
491  fixit(); /* use a macro for the 5 and 9 that appear here. */
492  /*We malloc the space here*/
493  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline =
494  (double *)MALLOC((unsigned long)(9)*sizeof(double));
495  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].SplineSecDer =
496  (double *)MALLOC((unsigned long)(9)*sizeof(double));
497 
498  /* always read at least 5 */
499  for( intsplinepts=0; intsplinepts<10; intsplinepts++ )
500  {
501  double temp = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
502 
503  if( lgEOL )
504  {
505  /* there should be 5 or 9 spline points. If we got EOL,
506  * insist that we were trying to read the 6th or 10th. */
507  if( (intsplinepts != 5) && (intsplinepts != 9) )
508  {
509  fprintf( ioQQQ, " There should have been 5 or 9 spline points. Sorry.\n" );
510  fprintf( ioQQQ, "string==%s==\n" ,chLine );
511  cdEXIT(EXIT_FAILURE);
512  }
513  else
514  break;
515  }
516  else
517  {
518  ASSERT( intsplinepts < 9 );
519  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[intsplinepts] = temp;
520  }
521  }
522 
523  ASSERT( (intsplinepts == 5) || (intsplinepts == 9) );
524 
525  /*The zeroth element contains the number of spline points*/
526  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].nSplinePts = intsplinepts;
527  /*Transition type*/
528  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].intTranType = intTranType;
529  /*gf value*/
530  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].gf = fGF;
531  /*Energy difference between two levels in Rydbergs*/
532  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].EnergyDiff = fEnergyDiff;
533  /*Scaling parameter C*/
534  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].ScalingParam = fScalingParam;
535 
536  /*Once the spline points have been filled,fill the second derivatives structure*/
537  /*Creating spline points array*/
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));
541 
542  if(intsplinepts == 5)
543  {
544  for(intxs=0;intxs<5;intxs++)
545  {
546  xs[intxs] = 0.25*intxs;
547  spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[intxs];
548  }
549  }
550  else if(intsplinepts == 9)
551  {
552  for(intxs=0;intxs<9;intxs++)
553  {
554  xs[intxs] = 0.125*intxs;
555  spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[intxs];
556  }
557  }
558  else
559  TotalInsanity();
560 
561  spline(xs, spl,intsplinepts,2e31,2e31,spl2);
562 
563  /*Filling the second derivative structure*/
564  for(i=0;i<intsplinepts;i++)
565  {
566  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].SplineSecDer[i] = spl2[i];
567  }
568 
569  free( xs );
570  free( spl );
571  free( spl2 );
572 
573  if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolEleColDATA ) == NULL )
574  {
575  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename);
576  cdEXIT(EXIT_FAILURE);
577  }
578  }
579 
580  /*****************************************/
581  /*** Print the AtmolCollSplines Values ***/
582  /*****************************************/
583 
584  if(DEBUGSTATE)
585  {
586  intNS = 0;
587  intCollIndex = 0;
588  ipHi = 1;
589  ipLo = 0;
590 
591  fprintf( ioQQQ, "Collisional Data: %li\t%li\t%f\t%f\t%f",
592  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].nSplinePts,
593  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].intTranType,
594  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].gf,
595  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].EnergyDiff,
596  AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].ScalingParam );
597  for( j=0; j< AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].nSplinePts; j++)
598  {
599  fprintf( ioQQQ, "\t%f",AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[j]);
600  }
601  fprintf( ioQQQ, "\n" );
602  }
603 
604  free( atmolStatesEnergy );
605  free( intNewIndex );
606 
607  fclose( atmolLevDATA );
608  fclose( atmolTraDATA );
609  fclose( atmolEleColDATA );
610 
611  return;
612 }

Generated for cloudy by doxygen 1.8.3.1