cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "lines_service.h" 00005 #include "physconst.h" 00006 #include "taulines.h" 00007 #include "trace.h" 00008 #include "string.h" 00009 #include "thirdparty.h" 00010 00011 #define DEBUGSTATE false 00012 00013 void atmdat_Chianti_readin( long intNS ) 00014 { 00015 DEBUG_ENTRY( "atmdat_Chianti_readin()" ); 00016 00017 int intCS,intCollIndex = -10000,intsplinepts,intTranType,intxs; 00018 /* type is set to 0 for non chianti and 1 for chianti*/ 00019 long int i,j,nMolLevs,intCollTran,intcollindex,ipLo,ipHi; 00020 FILE *atmolLevDATA , *atmolTraDATA=NULL,*atmolEleColDATA=NULL,*atmolProColDATA=NULL; 00021 realnum fstatwt,fenergyK,fenergyWN,fWLAng,fenergy,feinsteina; 00022 double fScalingParam,fEnergyDiff,fGF,*xs,*spl,*spl2; 00023 00024 //Energy Column Start 00025 //Statistical Weight Start 00026 //Lower Level Start 00027 //Upper Level Start 00028 //Einstein A Start 00029 //Wavelength of transition start 00030 long ECS,SWS,LLS,ULS,EAS,WOTS; 00031 00032 char chLine[FILENAME_PATH_LENGTH_2] , 00033 chEnFilename[FILENAME_PATH_LENGTH_2], 00034 chTraFilename[FILENAME_PATH_LENGTH_2], 00035 chEleColFilename[FILENAME_PATH_LENGTH_2], 00036 chProColFilename[FILENAME_PATH_LENGTH_2]; 00037 00038 char *chDesired; 00039 realnum *atmolStatesEnergy; 00040 long int *intNewIndex=NULL,*intOldIndex=NULL; 00041 int interror; 00042 bool lgProtonData=false,lgEneLevOrder; 00043 00044 /*This is to identify where the fields start */ 00045 /*In Chianti(1),the energy field in cm-1 starts at 40 and has a width 15*/ 00046 ECS = 71;/*We use the theoretical energy in cm-1*/ 00047 SWS = 38; 00048 LLS = 1; 00049 ULS = 6; 00050 EAS = 41; 00051 WOTS = 11; 00052 00053 /*For the CHIANTI DATABASE*/ 00054 /*First Opening the energy levels file*/ 00055 strcpy( chEnFilename , Species[intNS].chptrSpName ); 00056 strcat( chEnFilename , ".elvlc"); 00057 uncaps( chEnFilename ); 00058 if(DEBUGSTATE) 00059 printf("The data file is %s \n",chEnFilename); 00060 00061 /*Open the files*/ 00062 if( trace.lgTrace ) 00063 fprintf( ioQQQ," moldat_readin opening %s:",chEnFilename); 00064 00065 atmolLevDATA = open_data( chEnFilename, "r" ); 00066 00067 /*Second:Opening the transition probabilities file*/ 00068 strcpy( chTraFilename , Species[intNS].chptrSpName ); 00069 strcat( chTraFilename , ".wgfa"); 00070 uncaps( chTraFilename ); 00071 if(DEBUGSTATE) 00072 printf("The data file is %s \n",chTraFilename); 00073 00074 /*Open the files*/ 00075 if( trace.lgTrace ) 00076 fprintf( ioQQQ," moldat_readin opening %s:",chTraFilename); 00077 00078 atmolTraDATA = open_data( chTraFilename, "r" ); 00079 00080 /*Third: Opening the electron collision data*/ 00081 strcpy( chEleColFilename , Species[intNS].chptrSpName ); 00082 strcat( chEleColFilename , ".splups"); 00083 uncaps( chEleColFilename ); 00084 if(DEBUGSTATE) 00085 printf("The data file is %s \n",chEleColFilename); 00086 /*Open the files*/ 00087 if( trace.lgTrace ) 00088 fprintf( ioQQQ," moldat_readin opening %s:",chEleColFilename); 00089 00090 atmolEleColDATA = open_data( chEleColFilename, "r" ); 00091 00092 /*Fourth: Opening the proton collision data*/ 00093 strcpy( chProColFilename , Species[intNS].chptrSpName ); 00094 strcat( chProColFilename , ".psplups"); 00095 uncaps( chProColFilename ); 00096 if(DEBUGSTATE) 00097 printf("The data file is %s \n",chProColFilename); 00098 /*Open the files*/ 00099 if( trace.lgTrace ) 00100 fprintf( ioQQQ," moldat_readin opening %s:",chProColFilename); 00101 /*We will set a flag here to indicate if the proton collision strengths are available */ 00102 if( ( atmolProColDATA = fopen( chProColFilename , "r" ) ) != NULL ) 00103 { 00104 lgProtonData = true; 00105 fclose( atmolProColDATA ); 00106 atmolProColDATA = NULL; 00107 } 00108 else 00109 { 00110 lgProtonData = false; 00111 } 00112 00113 nMolLevs = 0; 00114 lgEneLevOrder = true; 00115 while( read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) != NULL ) 00116 { 00117 chDesired = strtok(chLine," \t"); 00118 intCS = atoi(chDesired); 00119 if (intCS == -1) 00120 { 00121 break; 00122 } 00123 else 00124 nMolLevs++; 00125 } 00126 00127 if(DEBUGSTATE) 00128 printf("The number of energy levels is %li",nMolLevs); 00129 00130 if( nMolLevs <= 0 ) 00131 { 00132 fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEnFilename ); 00133 fprintf( ioQQQ, "The file must be corrupted.\n" ); 00134 cdEXIT( EXIT_FAILURE ); 00135 } 00136 00137 Species[intNS].numLevels_max = nMolLevs; 00138 Species[intNS].numLevels_local = nMolLevs; 00139 /*Malloc the energy array to check if the energy levels are in order*/ 00140 atmolStatesEnergy = (realnum*)MALLOC((unsigned long)(nMolLevs)*sizeof(realnum)); 00141 /*malloc the States array*/ 00142 atmolStates[intNS] = (quantumState *)MALLOC((size_t)(nMolLevs)*sizeof(quantumState)); 00143 /*malloc the Transition array*/ 00144 atmolTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)nMolLevs); 00145 atmolTrans[intNS][0] = NULL; 00146 for( ipHi = 1; ipHi < nMolLevs; ipHi++) 00147 { 00148 atmolTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi); 00149 for(ipLo = 0; ipLo < ipHi; ipLo++) 00150 { 00151 atmolTrans[intNS][ipHi][ipLo].Lo = &atmolStates[intNS][ipLo]; 00152 atmolTrans[intNS][ipHi][ipLo].Hi = &atmolStates[intNS][ipHi]; 00153 atmolTrans[intNS][ipHi][ipLo].Emis = &DummyEmis; 00154 } 00155 } 00156 00157 for( intcollindex = 0; intcollindex<NUM_COLLIDERS; intcollindex++ ) 00158 { 00159 CollRatesArray[intNS][intcollindex] = NULL; 00160 } 00161 /*Allocating space just for the electron*/ 00162 CollRatesArray[intNS][0] = (double**)MALLOC((nMolLevs)*sizeof(double*)); 00163 for( ipHi = 0; ipHi<nMolLevs; ipHi++ ) 00164 { 00165 CollRatesArray[intNS][0][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double)); 00166 for( ipLo = 0; ipLo<nMolLevs; ipLo++) 00167 { 00168 CollRatesArray[intNS][0][ipHi][ipLo] = 0.0; 00169 } 00170 } 00171 /*Allocating space for the proton*/ 00172 if(lgProtonData) 00173 { 00174 CollRatesArray[intNS][1] = (double**)MALLOC((nMolLevs)*sizeof(double*)); 00175 for( ipHi = 0; ipHi<nMolLevs; ipHi++ ) 00176 { 00177 CollRatesArray[intNS][1][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double)); 00178 for( ipLo = 0; ipLo<nMolLevs; ipLo++) 00179 { 00180 CollRatesArray[intNS][1][ipHi][ipLo] = 0.0; 00181 } 00182 } 00183 } 00184 /*Rewind the energy levels files*/ 00185 if( fseek( atmolLevDATA , 0 , SEEK_SET ) != 0 ) 00186 { 00187 fprintf( ioQQQ, " moldat_readin could not rewind %s.\n",chEnFilename); 00188 cdEXIT(EXIT_FAILURE); 00189 } 00190 00191 00192 /*Reading in the energy and checking that they are in order*/ 00193 for( i=0; i<nMolLevs; i++ ) 00194 { 00195 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00196 { 00197 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename); 00198 cdEXIT(EXIT_FAILURE); 00199 } 00200 fenergy = (realnum) atof(&chLine[ECS]); 00201 atmolStatesEnergy[i] = fenergy; 00202 00203 /*To check if energy levels are in order*/ 00204 /*If the energy levels are not in order a flag is set*/ 00205 if(i>0) 00206 { 00207 if (atmolStatesEnergy[i] < atmolStatesEnergy[i-1]) 00208 { 00209 lgEneLevOrder = false; 00210 } 00211 } 00212 } 00213 00214 if(DEBUGSTATE) 00215 { 00216 for(i=0;i<nMolLevs;i++) 00217 { 00218 printf("The %ld value of the unsorted energy array is %f \n",i,atmolStatesEnergy[i]); 00219 } 00220 } 00221 00222 intNewIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int)); 00223 00224 if(lgEneLevOrder == false) 00225 { 00226 /*If the energy levels are not in order and the flag is set(lgEneLevOrder=FALSE) 00227 then we sort the energy levels to find the relation matrix between the old and new indices*/ 00228 intOldIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int)); 00229 /*We now do a modified quick sort*/ 00230 spsort(atmolStatesEnergy,nMolLevs,intOldIndex,2,&interror); 00231 /*intNewIndex has the key to map*/ 00232 for( i=0; i<nMolLevs; i++ ) 00233 { 00234 ASSERT( intOldIndex[i] < nMolLevs ); 00235 intNewIndex[intOldIndex[i]] = i; 00236 } 00237 if(DEBUGSTATE) 00238 { 00239 for(i=0;i<nMolLevs;i++) 00240 { 00241 printf("The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]); 00242 } 00243 for( i=0; i<nMolLevs; i++ ) 00244 { 00245 printf("The %ld value of the sorted energy array is %f \n",i,atmolStatesEnergy[i]); 00246 } 00247 } 00248 00249 free( intOldIndex ); 00250 } 00251 else 00252 { 00253 /* if energies already in correct order, new index is same as original. */ 00254 for( i=0; i<nMolLevs; i++ ) 00255 { 00256 intNewIndex[i] = i; 00257 } 00258 } 00259 00260 /* insist that there the intNewIndex values are unique */ 00261 for( i=0; i<nMolLevs; i++ ) 00262 { 00263 for( j=i+1; j<nMolLevs; j++ ) 00264 { 00265 ASSERT( intNewIndex[i] != intNewIndex[j] ); 00266 } 00267 } 00268 00269 00270 /*After we read in the energies we rewind the energy levels file */ 00271 if( fseek( atmolLevDATA , 0 , SEEK_SET ) != 0 ) 00272 { 00273 fprintf( ioQQQ, " moldat_readin could not rewind %s.\n",chEnFilename); 00274 cdEXIT(EXIT_FAILURE); 00275 } 00276 00277 for( i=0; i<nMolLevs; i++ ) 00278 { 00279 long j = intNewIndex[i]; 00280 00281 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00282 { 00283 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename); 00284 cdEXIT(EXIT_FAILURE); 00285 } 00286 /*information needed for label*/ 00287 strcpy(atmolStates[intNS][j].chLabel, " "); 00288 strncpy(atmolStates[intNS][j].chLabel,Species[intNS].chptrSpName, 4); 00289 atmolStates[intNS][j].chLabel[4] = '\0'; 00290 00291 fstatwt = (realnum)atof(&chLine[SWS]); 00292 if(fstatwt <= 0.) 00293 { 00294 fprintf( ioQQQ, " A positive non zero value is expected for the statistical weight .\n"); 00295 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename); 00296 cdEXIT(EXIT_FAILURE); 00297 00298 } 00299 atmolStates[intNS][j].g = fstatwt; 00300 fenergy = (realnum) atof(&chLine[ECS]); 00301 atmolStates[intNS][j].energy = fenergy; 00302 00303 if(DEBUGSTATE) 00304 { 00305 fprintf( ioQQQ, "The %ld value of g after populating is %f \n", 00306 j,atmolStates[intNS][j].g); 00307 fprintf( ioQQQ, "The %ld value of energy after populating is %f \n", 00308 j,atmolStates[intNS][j].energy); 00309 } 00310 } 00311 00312 00313 /* fill in all transition energies, can later overwrite for specific radiative transitions */ 00314 for( i=1; i<nMolLevs; i++ ) 00315 { 00316 for( j=0; j<i; j++ ) 00317 { 00318 fenergyWN = MAX2( (realnum)1E-10, atmolStates[intNS][i].energy - atmolStates[intNS][j].energy ); 00319 fenergyK = (realnum)(fenergyWN*T1CM); 00320 00321 atmolTrans[intNS][i][j].EnergyK = fenergyK; 00322 atmolTrans[intNS][i][j].EnergyWN = fenergyWN; 00323 atmolTrans[intNS][i][j].EnergyErg = (realnum)ERG1CM *fenergyWN; 00324 atmolTrans[intNS][i][j].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN)); 00325 } 00326 } 00327 00328 /************************************************************************/ 00329 /*** Read in the level numbers, Einstein A and transition wavelength ***/ 00330 /************************************************************************/ 00331 if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolTraDATA ) == NULL ) 00332 { 00333 fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename); 00334 cdEXIT(EXIT_FAILURE); 00335 } 00336 00337 while(atoi(chLine)!= -1) 00338 { 00339 ipLo = atoi(&chLine[LLS])-1; 00340 ipHi = atoi(&chLine[ULS])-1; 00341 00342 /* translate to properly sorted indices */ 00343 ipLo = intNewIndex[ipLo]; 00344 ipHi = intNewIndex[ipHi]; 00345 00346 ASSERT( ipHi > ipLo ); 00347 00348 fWLAng = (realnum)atof(&chLine[WOTS]); 00349 00350 /* \todo 2 Chianti labels the H 1 2-photon transition as z wavelength of zero. 00351 * Should we just ignore all of the wavelengths in this file and use the 00352 * difference of level energies instead. */ 00353 if( fWLAng == 0. ) 00354 { 00355 fWLAng = (realnum)1e8/ 00356 (atmolStates[intNS][ipHi].energy - atmolStates[intNS][ipLo].energy); 00357 } 00358 00359 feinsteina = (realnum)atof(&chLine[EAS]); 00360 atmolTrans[intNS][ipHi][ipLo].Emis 00361 = AddLine2Stack(feinsteina, &atmolTrans[intNS][ipHi][ipLo]); 00362 atmolTrans[intNS][ipHi][ipLo].WLAng = fWLAng; 00363 fenergyWN = (realnum)(1e+8/fWLAng); 00364 fenergyK = (realnum)(fenergyWN*T1CM); 00365 00366 atmolTrans[intNS][ipHi][ipLo].EnergyK = fenergyK; 00367 atmolTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN; 00368 atmolTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN; 00369 atmolTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN)); 00370 00371 ASSERT( !isnan( atmolTrans[intNS][ipHi][ipLo].EnergyK ) ); 00372 00373 if(DEBUGSTATE) 00374 { 00375 fprintf( ioQQQ, "The lower level is %ld \n",ipLo); 00376 fprintf( ioQQQ, "The upper level is %ld \n",ipHi); 00377 fprintf( ioQQQ, "The Einstein A is %f \n",atmolTrans[intNS][ipHi][ipLo].Emis->Aul); 00378 fprintf( ioQQQ, "The wavelength of the transition is %f \n",atmolTrans[intNS][ipHi][ipLo].WLAng ); 00379 } 00380 00381 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolTraDATA ) == NULL ) 00382 { 00383 fprintf( ioQQQ, " Data is expected on this line.\n"); 00384 fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename); 00385 cdEXIT(EXIT_FAILURE); 00386 } 00387 } 00388 00389 /* Malloc space for splines */ 00390 00391 AtmolCollSplines[intNS] = (CollSplinesArray***)MALLOC(nMolLevs *sizeof(CollSplinesArray**)); 00392 for( i=0; i<nMolLevs; i++ ) 00393 { 00394 AtmolCollSplines[intNS][i] = (CollSplinesArray **)MALLOC((unsigned long)(nMolLevs)*sizeof(CollSplinesArray *)); 00395 for( j=0; j<nMolLevs; j++ ) 00396 { 00397 AtmolCollSplines[intNS][i][j] = 00398 (CollSplinesArray *)MALLOC((unsigned long)(NUM_COLLIDERS)*sizeof(CollSplinesArray )); 00399 00400 for( long k=0; k<NUM_COLLIDERS; k++ ) 00401 { 00402 /* initialize all spline variables */ 00403 AtmolCollSplines[intNS][i][j][k].collspline = NULL; 00404 AtmolCollSplines[intNS][i][j][k].SplineSecDer = NULL; 00405 AtmolCollSplines[intNS][i][j][k].nSplinePts = -1; 00406 AtmolCollSplines[intNS][i][j][k].intTranType = -1; 00407 AtmolCollSplines[intNS][i][j][k].gf = BIGDOUBLE; 00408 AtmolCollSplines[intNS][i][j][k].EnergyDiff = BIGDOUBLE; 00409 AtmolCollSplines[intNS][i][j][k].ScalingParam = BIGDOUBLE; 00410 } 00411 } 00412 } 00413 00414 fixit(); /* This is where the loop of colliders should begin */ 00415 00416 /*********************************************/ 00417 /*** Read in the electron collisional data ***/ 00418 /*********************************************/ 00419 intCollIndex = 0; 00420 intCollTran = 0; 00421 if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolEleColDATA ) == NULL ) 00422 { 00423 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename); 00424 cdEXIT(EXIT_FAILURE); 00425 } 00426 /*We need to find the number of collisional transitions */ 00427 while(atoi(chLine)!= -1) 00428 { 00429 intCollTran++; 00430 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolEleColDATA ) == NULL ) 00431 { 00432 fprintf( ioQQQ, " Data is expected on this line.\n"); 00433 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename); 00434 cdEXIT(EXIT_FAILURE); 00435 00436 } 00437 } 00438 00439 if(DEBUGSTATE) 00440 printf("\n The number of collisional transitions is %li",intCollTran); 00441 00442 /*We rewind the file*/ 00443 if( fseek( atmolEleColDATA , 0 , SEEK_SET ) != 0 ) 00444 { 00445 fprintf( ioQQQ, " moldat_readin could not rewind %s.\n",chEleColFilename); 00446 cdEXIT(EXIT_FAILURE); 00447 } 00448 /*Dummy string used for convenience*/ 00449 strcpy(chLine,"A"); 00450 00451 if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolEleColDATA ) == NULL ) 00452 { 00453 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename); 00454 cdEXIT(EXIT_FAILURE); 00455 } 00456 00457 /*If the file exists,the first line cannot be "-1"*/ 00458 if(atoi(chLine)== -1) 00459 { 00460 fprintf( ioQQQ, " If the file %s exists,the first line cannot be -1 .\n",chEleColFilename); 00461 cdEXIT(EXIT_FAILURE); 00462 } 00463 00464 while(atoi(chLine)!= -1) 00465 { 00466 i = 1; 00467 bool lgEOL; 00468 00469 (void)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); // intAtomicNo 00470 (void)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); // intIonStage 00471 /* level indices */ 00472 ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1; 00473 ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1; 00474 intTranType = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00475 /*gf*/ 00476 fGF = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00477 /*Energy Difference*/ 00478 fEnergyDiff = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00479 /*Scaling Parameter*/ 00480 fScalingParam = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00481 00482 /* translate to properly sorted indices */ 00483 ipLo = intNewIndex[ipLo]; 00484 ipHi = intNewIndex[ipHi]; 00485 00486 ASSERT( ipLo < nMolLevs ); 00487 ASSERT( ipHi < nMolLevs ); 00488 ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline == NULL ); 00489 ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].SplineSecDer == NULL ); 00490 00491 fixit(); /* use a macro for the 5 and 9 that appear here. */ 00492 /*We malloc the space here*/ 00493 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline = 00494 (double *)MALLOC((unsigned long)(9)*sizeof(double)); 00495 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].SplineSecDer = 00496 (double *)MALLOC((unsigned long)(9)*sizeof(double)); 00497 00498 /* always read at least 5 */ 00499 for( intsplinepts=0; intsplinepts<10; intsplinepts++ ) 00500 { 00501 double temp = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00502 00503 if( lgEOL ) 00504 { 00505 /* there should be 5 or 9 spline points. If we got EOL, 00506 * insist that we were trying to read the 6th or 10th. */ 00507 if( (intsplinepts != 5) && (intsplinepts != 9) ) 00508 { 00509 fprintf( ioQQQ, " There should have been 5 or 9 spline points. Sorry.\n" ); 00510 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00511 cdEXIT(EXIT_FAILURE); 00512 } 00513 else 00514 break; 00515 } 00516 else 00517 { 00518 ASSERT( intsplinepts < 9 ); 00519 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[intsplinepts] = temp; 00520 } 00521 } 00522 00523 ASSERT( (intsplinepts == 5) || (intsplinepts == 9) ); 00524 00525 /*The zeroth element contains the number of spline points*/ 00526 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].nSplinePts = intsplinepts; 00527 /*Transition type*/ 00528 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].intTranType = intTranType; 00529 /*gf value*/ 00530 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].gf = fGF; 00531 /*Energy difference between two levels in Rydbergs*/ 00532 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].EnergyDiff = fEnergyDiff; 00533 /*Scaling parameter C*/ 00534 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].ScalingParam = fScalingParam; 00535 00536 /*Once the spline points have been filled,fill the second derivatives structure*/ 00537 /*Creating spline points array*/ 00538 xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00539 spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00540 spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double)); 00541 00542 if(intsplinepts == 5) 00543 { 00544 for(intxs=0;intxs<5;intxs++) 00545 { 00546 xs[intxs] = 0.25*intxs; 00547 spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[intxs]; 00548 } 00549 } 00550 else if(intsplinepts == 9) 00551 { 00552 for(intxs=0;intxs<9;intxs++) 00553 { 00554 xs[intxs] = 0.125*intxs; 00555 spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[intxs]; 00556 } 00557 } 00558 else 00559 TotalInsanity(); 00560 00561 spline(xs, spl,intsplinepts,2e31,2e31,spl2); 00562 00563 /*Filling the second derivative structure*/ 00564 for(i=0;i<intsplinepts;i++) 00565 { 00566 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].SplineSecDer[i] = spl2[i]; 00567 } 00568 00569 free( xs ); 00570 free( spl ); 00571 free( spl2 ); 00572 00573 if(read_whole_line( chLine , (int)sizeof(chLine) ,atmolEleColDATA ) == NULL ) 00574 { 00575 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEleColFilename); 00576 cdEXIT(EXIT_FAILURE); 00577 } 00578 } 00579 00580 /*****************************************/ 00581 /*** Print the AtmolCollSplines Values ***/ 00582 /*****************************************/ 00583 00584 if(DEBUGSTATE) 00585 { 00586 intNS = 0; 00587 intCollIndex = 0; 00588 ipHi = 1; 00589 ipLo = 0; 00590 00591 fprintf( ioQQQ, "Collisional Data: %li\t%li\t%f\t%f\t%f", 00592 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].nSplinePts, 00593 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].intTranType, 00594 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].gf, 00595 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].EnergyDiff, 00596 AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].ScalingParam ); 00597 for( j=0; j< AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].nSplinePts; j++) 00598 { 00599 fprintf( ioQQQ, "\t%f",AtmolCollSplines[intNS][ipHi][ipLo][intCollIndex].collspline[j]); 00600 } 00601 fprintf( ioQQQ, "\n" ); 00602 } 00603 00604 free( atmolStatesEnergy ); 00605 free( intNewIndex ); 00606 00607 fclose( atmolLevDATA ); 00608 fclose( atmolTraDATA ); 00609 fclose( atmolEleColDATA ); 00610 00611 return; 00612 }