root/branches/newmole/source/atmdat_lamda.cpp

Revision 2593, 14.8 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
13/*Separate Routine to read in the molecules*/
14void atmdat_lamda_readin( long intNS )
15{
16        DEBUG_ENTRY( "atmdat_lamda_readin()" );
17
18        int ngMolLevs = -10000,intCollIndex = -10000,intLColliderIndex = -10000;
19        /* type is set to 0 for non chianti and 1 for chianti*/
20        long int i,j,nMolLevs,intlnct,intrtct,intgrtct,intCollPart,
21                intDCollPart,intCollTran,nCollTrans,intCollTemp,intcollindex,
22                ipLo,ipHi;
23        /*intrtct refers to radiative transitions count*/
24        FILE *atmolLevDATA;
25        realnum  fstatwt,fenergyK,fenergyWN,fenergy,feinsteina;
26
27        char chLine[FILENAME_PATH_LENGTH_2] ,
28                chEFilename[FILENAME_PATH_LENGTH_2],
29                *chColltemp,*chCollName;
30        char chDirectory[10];
31
32
33        ASSERT( lgSpeciesMolecule[intNS] );
34
35        /*Load the species name in the Species array structure*/
36        if(DEBUGSTATE)
37                printf("The name of the %li species is %s \n",intNS+1,Species[intNS].chptrSpName);
38
39#       ifdef _WIN32
40        strcpy( chDirectory, "lamda\\" );
41#       else
42        strcpy( chDirectory, "lamda/" );
43#       endif
44
45        strcpy( chEFilename , chDirectory );
46        strcat( chEFilename , Species[intNS].chptrSpName );                     
47        strcat( chEFilename , ".dat");
48        uncaps( chEFilename );
49
50        /*Open the files*/
51        if( trace.lgTrace )
52                fprintf( ioQQQ," moldat_readin opening %s:",chEFilename);
53
54        atmolLevDATA = open_data( chEFilename, "r" );
55
56        nMolLevs = 0;
57        ngMolLevs = 0;
58        intrtct = 0;
59        intgrtct = 0;
60        intlnct = 0;
61        ipHi = 0;
62        ipLo = 0;
63        j = 0;
64        while(intlnct < 3)
65        {
66                intlnct++;
67                if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
68                {
69                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
70                        cdEXIT(EXIT_FAILURE);
71                }
72        }
73        /*Extracting out the molecular weight*/
74        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
75        {
76                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
77                cdEXIT(EXIT_FAILURE);
78        }
79        Species[intNS].fmolweight = (realnum)atof(chLine);
80
81        /*Discard this line*/
82        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
83        {
84                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
85                cdEXIT(EXIT_FAILURE);
86        }
87        /*Reading in the number of energy levels*/
88        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
89        {
90                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
91                cdEXIT(EXIT_FAILURE);
92        }
93        ngMolLevs = atoi(chLine);
94
95        if(ngMolLevs <= 0)
96        {
97                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
98                cdEXIT(EXIT_FAILURE);
99        }
100        Species[intNS].numLevels_max = ngMolLevs;
101        Species[intNS].numLevels_local = ngMolLevs;
102        if(DEBUGSTATE)
103        {
104                printf("The number of energy levels is %li \n",Species[intNS].numLevels_max);
105        }
106        /*the above refers to the number of energy levels specified in the data file*/
107        /*malloc the States array*/
108        atmolStates[intNS] = (quantumState *)MALLOC((size_t)(ngMolLevs)*sizeof(quantumState));
109        /*malloc the Transition array*/
110        atmolTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)ngMolLevs);
111        atmolTrans[intNS][0] = NULL;
112        for( ipHi = 1; ipHi < ngMolLevs; ipHi++)
113        {
114                atmolTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi);
115                for(ipLo = 0; ipLo < ipHi; ipLo++)
116                {
117                        atmolTrans[intNS][ipHi][ipLo].Lo = &atmolStates[intNS][ipLo];
118                        atmolTrans[intNS][ipHi][ipLo].Hi = &atmolStates[intNS][ipHi];
119                        atmolTrans[intNS][ipHi][ipLo].Emis = &DummyEmis;
120                }
121        }
122
123        for( intcollindex = 0; intcollindex <NUM_COLLIDERS; intcollindex++ )
124        {
125                CollRatesArray[intNS][intcollindex] = NULL;
126        }
127
128
129        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
130        {
131                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
132                cdEXIT(EXIT_FAILURE);
133        }
134
135        for( nMolLevs=0; nMolLevs<ngMolLevs; nMolLevs++)
136        {       
137                if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
138                {
139                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
140                        cdEXIT(EXIT_FAILURE);
141                }
142                /*information needed for label*/
143                strcpy(atmolStates[intNS][nMolLevs].chLabel, "    ");
144                strncpy(atmolStates[intNS][nMolLevs].chLabel,Species[intNS].chptrSpName, 4);
145                atmolStates[intNS][nMolLevs].chLabel[4] = '\0';
146
147                long i = 1;
148                bool lgEOL;
149                long index;
150
151                index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
152                fenergy = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
153                fstatwt = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
154
155                ASSERT( index == nMolLevs + 1 );
156
157                atmolStates[intNS][nMolLevs].energy = fenergy;
158                atmolStates[intNS][nMolLevs].g = fstatwt;
159
160                if (nMolLevs > 0)
161                {
162                        if (atmolStates[intNS][nMolLevs].energy < atmolStates[intNS][nMolLevs-1].energy)
163                        {
164                                fprintf( ioQQQ, " The energy levels are not in order");
165                                cdEXIT(EXIT_FAILURE);
166                        }
167                }
168                if(DEBUGSTATE)
169                {
170                        printf("The converted energy is %f \n",atmolStates[intNS][nMolLevs].energy);
171                        printf("The value of g is %f \n",atmolStates[intNS][nMolLevs].g);
172                }
173        }
174
175        /* fill in all transition energies, can later overwrite for specific radiative transitions */
176        for( i=1; i<ngMolLevs; i++ )
177        {
178                for( j=0; j<i; j++ )
179                {
180                        fenergyWN = atmolStates[intNS][i].energy - atmolStates[intNS][j].energy;
181                        fenergyK = (realnum)(fenergyWN*T1CM);
182
183                        atmolTrans[intNS][i][j].EnergyK = fenergyK;
184                        atmolTrans[intNS][i][j].EnergyWN = fenergyWN;
185                        atmolTrans[intNS][i][j].EnergyErg = (realnum)ERG1CM *fenergyWN;
186
187                        /* there are OH hyperfine levels where i+1 and i have exactly
188                         * the same energy.  The refractive index routine will FPE with
189                         * an energy of zero - so we do this test */
190                        if( fenergyWN>SMALLFLOAT )
191                                atmolTrans[intNS][i][j].WLAng = (realnum)(1e+8f/fenergyWN/
192                                        RefIndex(fenergyWN));
193                        else
194                                atmolTrans[intNS][i][j].WLAng = 1e30f;
195                }
196        }
197
198        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
199        {
200                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
201                cdEXIT(EXIT_FAILURE);
202        }
203        if(chLine[0]!='!')
204        {
205                fprintf( ioQQQ, " The number of energy levels in file %s is not correct.\n",chEFilename);
206                cdEXIT(EXIT_FAILURE);
207        }
208        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
209        {
210                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
211                cdEXIT(EXIT_FAILURE);
212        }
213        intgrtct = atoi(chLine);
214        /*The above gives the number of radiative transitions*/
215        if(intgrtct <= 0)
216        {
217                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
218                cdEXIT(EXIT_FAILURE);
219        }
220        if(DEBUGSTATE)
221        {
222                printf("The number of radiative transitions is %li \n",intgrtct);
223        }
224        if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
225        {
226                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
227                cdEXIT(EXIT_FAILURE);
228        }
229
230        for( intrtct=0; intrtct<intgrtct; intrtct++)
231        {
232                if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
233                {
234                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
235                        cdEXIT(EXIT_FAILURE);
236                }
237
238                long i = 1;
239                bool lgEOL;
240                long index;
241
242                index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
243                ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
244                ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
245                feinsteina = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
246                /* don't need the energy in GHz, so throw it away. */
247                FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
248                fenergyK = (realnum)((atmolStates[intNS][ipHi].energy -atmolStates[intNS][ipLo].energy)*T1CM);
249                ASSERT( index == intrtct + 1 );
250
251                atmolTrans[intNS][ipHi][ipLo].Emis = AddLine2Stack(feinsteina, &atmolTrans[intNS][ipHi][ipLo]);
252                atmolTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
253                fenergyWN = (realnum)((fenergyK)/T1CM);
254                atmolTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
255                atmolTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
256                atmolTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
257
258                ASSERT( !isnan( atmolTrans[intNS][ipHi][ipLo].EnergyK ) );
259
260                if(DEBUGSTATE)
261                {
262                        printf("The upper level is %ld \n",ipHi+1);
263                        printf("The lower level is %ld \n",ipLo+1);
264                        printf("The Einstein A  is %E \n",atmolTrans[intNS][ipHi][ipLo].Emis->Aul);
265                        printf("The Energy of the transition is %E \n",atmolTrans[intNS][ipHi][ipLo].EnergyK);
266                }
267        }
268
269        if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
270        {
271                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
272                cdEXIT(EXIT_FAILURE);
273        }
274
275        if(chLine[0]!='!')
276        {
277                fprintf( ioQQQ, " The number of radiative transitions in file %s is not correct.\n",chEFilename);
278                cdEXIT(EXIT_FAILURE);
279        }
280
281        /*Getting the number of collisional partners*/
282        if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
283        {
284                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
285                cdEXIT(EXIT_FAILURE);
286        }
287        else
288        {
289                intCollPart = atoi(chLine);
290        }
291        /*Checking the number of collisional partners does not exceed 9*/
292        if(intCollPart > NUM_COLLIDERS-1)
293        {
294                fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename );
295                cdEXIT(EXIT_FAILURE);
296        }
297        /*Creating the duplicate of the number of collisional partners which is reduced*/
298        intDCollPart = intCollPart;
299        while(intDCollPart > 0)
300        {
301                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
302                {
303                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
304                        cdEXIT(EXIT_FAILURE);
305                }
306
307                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
308                {
309                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
310                        cdEXIT(EXIT_FAILURE);
311                }
312                /*Extract out the name of the collider*/
313                /*The following are the rules expected in the datafiles to extract the names of the collider*/
314                /*The line which displays the species and the collider starts with a number*/
315                /*This refers to the collider in the Leiden database*/
316                /*In the Leiden database 1 referes to H2,2 to para-H2,3 to ortho-H2
317                4 to electrons,5 to H and 6 to He*/
318                chCollName = strtok(chLine," ");
319                /*Leiden Collider Index*/
320                intLColliderIndex = atoi(chCollName);
321                /*In Cloudy,We assign the following indices for the colliders*/
322                /*electron=0,proton=1,atomic hydrogen=2,He=3,He+=4,He++=5,oH2=6,pH2=7,H2=8*/
323
324                if(intLColliderIndex == 1)
325                {
326                        intCollIndex = 8;
327                }
328                else if(intLColliderIndex == 2)
329                {
330                        intCollIndex = 7;
331                }
332                else if(intLColliderIndex == 3)
333                {
334                        intCollIndex = 6;
335                }
336                else if(intLColliderIndex == 4)
337                {
338                        intCollIndex = 0;
339                }
340                else if(intLColliderIndex == 5)
341                {
342                        intCollIndex = 2;
343                }
344                else if(intLColliderIndex == 6)
345                {
346                        intCollIndex = 3;
347                }
348                else
349                {
350                        TotalInsanity();
351                }
352                /*This is where we allocate memory if the collider exists*/
353                /*Needed to take care of the he collisions*/
354                CollRatesArray[intNS][intCollIndex] = (double**)MALLOC((unsigned long)(ngMolLevs) * sizeof(double*));
355                for( ipHi = 0; ipHi<ngMolLevs; ipHi++ )
356                {
357                        CollRatesArray[intNS][intCollIndex][ipHi] = (double*)MALLOC((unsigned long)(ngMolLevs) * sizeof(double));
358                        for( ipLo = 0; ipLo<ngMolLevs; ipLo++ )
359                        {
360                                CollRatesArray[intNS][intCollIndex][ipHi][ipLo] = 0.0;
361                        }
362                }
363                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
364                {
365                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
366                        cdEXIT(EXIT_FAILURE);
367                }
368                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
369                {
370                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
371                        cdEXIT(EXIT_FAILURE);
372                }
373                /*Number of coll trans*/
374                intCollTran = atoi(chLine);
375
376                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
377                {
378                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
379                        cdEXIT(EXIT_FAILURE);
380                }
381                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
382                {
383                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
384                        cdEXIT(EXIT_FAILURE);
385                }
386                /*Number of coll temps*/
387                intCollTemp = atoi(chLine);
388                /*Storing the number of collisional temperatures*/
389                AtmolCollRateCoeff[intNS][intCollIndex].ntemps = intCollTemp;
390                /*Mallocing*/
391                AtmolCollRateCoeff[intNS][intCollIndex].temps =
392                        (double *)MALLOC((unsigned long)intCollTemp*sizeof(double));
393                AtmolCollRateCoeff[intNS][intCollIndex].collrates =
394                        (double***)MALLOC((unsigned long)(ngMolLevs)*sizeof(double**));
395                for( i=0; i<ngMolLevs; i++ )
396                {
397                        AtmolCollRateCoeff[intNS][intCollIndex].collrates[i] =
398                                (double **)MALLOC((unsigned long)(ngMolLevs)*sizeof(double*));
399                        for( j=0; j<ngMolLevs; j++ )
400                        {
401                                AtmolCollRateCoeff[intNS][intCollIndex].collrates[i][j] =
402                                        (double *)MALLOC((unsigned long)(intCollTemp)*sizeof(double));
403                        }
404                }
405
406                /*Discard the header line*/
407                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
408                {
409                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
410                        cdEXIT(EXIT_FAILURE);
411                }
412                /*Getting the collisional Temperatures*/
413                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
414                {
415                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
416                        cdEXIT(EXIT_FAILURE);
417                }
418                /*Filling the collisional temps array*/
419                chColltemp =strtok(chLine," ");
420                AtmolCollRateCoeff[intNS][intCollIndex].temps[0] =(realnum) atof(chColltemp);
421                for( i=1; i<intCollTemp; i++ )
422                {
423                        chColltemp =strtok(NULL," ");
424                        AtmolCollRateCoeff[intNS][intCollIndex].temps[i] = atof(chColltemp);
425                }
426                /*Discard the header line*/
427                if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
428                {
429                        fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
430                        cdEXIT(EXIT_FAILURE);
431                }
432                /*Getting all the collisional transitions data*/
433                for( nCollTrans=0; nCollTrans<intCollTran; nCollTrans++ )
434                {
435                        if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
436                        {
437                                fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
438                                cdEXIT(EXIT_FAILURE);
439                        }
440
441                        long i = 1;
442                        bool lgEOL;
443                        long index;
444
445                        index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
446                        ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
447                        ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
448
449                        /* Indices between the very highest levels seem to be reversed */
450                        if( ipHi < ipLo )
451                        {
452                                ASSERT( ipLo == ngMolLevs - 1);
453                                long temp = ipHi;
454                                ipHi = ipLo;
455                                ipLo = temp;
456                        }
457
458                        for( long j=0; j<intCollTemp; j++ )
459                        {
460                                AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][j] =
461                                        (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
462                        }
463
464                        ASSERT( index == nCollTrans + 1 );
465
466                        if(DEBUGSTATE)
467                        {
468                                printf("The values of up and lo are %ld & %ld \n",ipHi,ipLo);
469                                printf("The collision rates are");
470                                for(i=0;i<intCollTemp;i++)
471                                {
472                                        printf("\n %e",AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][i]);
473                                }
474                                printf("\n");
475                        }
476                }
477
478                intDCollPart = intDCollPart -1;
479        }
480
481        return;
482}
Note: See TracBrowser for help on using the browser.