root/branches/newmole/source/atmdat_chianti.cpp

Revision 2593, 19.4 kB (checked in by rporter, 13 days ago)

branches/newmole - bring in trunk r2581:2592

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